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ABSTRACT 


This report relates the shape of the retarding potential curve produced by gridded 
electrostatic analyzers to the floating potential, the kinetic temperature, and the half- 
angle subtended in velocity space by the particles. Special emphasis is placed on the 
case in which the particles are Maxwellianiz ed, but restricted within an ^ ^escape cone” 
in velocity space of half-angle 6 ^, The cases of pulsed and steady- state supply of parti- 
cles are both considered. A computer program is included which produces an iterated 
best fit of experimental integrated energy spectra to these analjrtical expressions, and 
provides the effective values of the kinetic temperature, floating potential, and escape 
cone angle. A series of experimental runs is analyzed with the computer program as an 
illustrative example, and the possible sources of error are discussed. 
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ANALYSIS OF INTEGRATED CHARGED PARTICLE ENERGY SPECTRA 
FROM GRIDDED ELECTROSTATIC ANALYZERS 
by J. Reece Roth and Marion Clark 
Lewis Research Center 

SUMMARY 

This report presents a systematic means of evaluating the retarding potential curves 
obtained from gridded electrostatic analyzers. Special emphasis is placed on the case in 
which the particles have a Maxwellian velocity distribution along a radius in velocity 
space, but are confined within an ’’escape cone” in velocity space of half- angle 0^. Such 
velocity distributions may result when particles are lost from ’’magnetic mirror” con- 
finement geometries and also when an isotropic and Maxwellian group of particles passes 
through a series of slits or apertures. 

Analytical expressions are derived for the expected integrated energy spectra as a 
function of the half- angle 0^. The two distinct cases of a pulsed and a steady- state sup- 
ply of particles are considered. The pulsed case would be appropriate to the analysis of 
a laser- generated plasma, to certain plasmas generated by rapid magnetic compression, 
or, under certain conditions, to the experiment discussed herein. The steady- state case 
is appropriate to some of the experiments performed to illustrate the data analysis tech- 
nique and to other situations of interest in controlled fusion research. 

A computer program is included which obtains an iterated best fit of experimental 
energy spectra to these analytical expressions. This computer program converges on 
values of the kinetic temperature, floating potential, and escape cone angle which best fit 
the experimental data. A series of illustrative experimental runs is reduced with the 
computer program, and the possible sources of error are discussed. By using the com- 
puter program to compare the analytical expressions derived herein with experimental 
retarding potential curves, significantly more information may reliably be obtained about 
the velocity distribution than is possible with conventional methods of data reduction, 
which generally yield only the kinetic temperature of the charged particles. 




INTRODUCTION 


Retarding field electrostatic analyzers have been widely used to determine the distri- 
bution of momentum in a flux of charged particles in the fields of plasma physics (refs. 1 
to 5), solid-state physics (refs. 6 and 7), and in other applications. The various types 
of retarding field energy analyzers have been discussed by Simpson (ref. 8), who points 
out their individual advantages and difficulties. The specific type of retarding field en- 
ergy analyzer used in the experimental portion of this investigation is the cylindrically 
symmetric gridded electrostatic analyzer illustrated in figure 1. In use, a flux of 
charged particles approaches the analyzer from the right. After passing through a 
grounded aperture plate, which establishes a zero reference for potential, the particles 
impinge on grid 3, which may be biased to discriminate between and repel all charges of 
either positive or negative sign. The particles then approach grid 2, which is biased to 
repel all particles whose momentum in the direction of the analyzer axis is below the 
value corresponding to the potential on this grid. The particles with momentum above 
this value pass through grid 2, and then pass through the negatively biased grid 1, whose 
function is to prevent secondary electrons from leaving the collector. Varying the bias 
voltage V on grid 2 from zero to a value much higher than that corresponding to the par- 
ticle kinetic temperature can make the current to the collector produce a retarding poten- 


/-Grid 2: 35 mesh 
/ wire fabric 



Figure 1. - Schematic drawing of gridded electrostatic analyzer used In 
experimental investigation. 
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tial curve of V as a function of collector current. This curve represents the integrated 
distribution of momentum of particles along the collector axis. Since the analyzer has 
been assumed to be one dimensional, it will provide information only about the momentum 
along the analyzer axis, and not necessarily about the total energy of the particles. 

hi order to establish a one-to-one relation between the particle kinetic temperature 
and the axial momentum, effects which alter the axial momentum of the particles after 
they enter the analyzer must be eliminated or accounted for. Much of the existing litera- 
ture has been devoted to the analysis of such effects. The focusing effect of grids and 
apertures on the axial momentum of the particles has been considered by Orlinov (ref. 7), 
Simpson (ref. 8), Simpson and Marton (ref. 9), and Mendlowitz (ref. 10). Space-charge 
and sheath effects can also alter the axial momentum of the particles if the plasma den- 
sity between the grids is sufficiently hi^. Such effects were considered by Mason (ref. 4) 
and by Fleischmann, et al. (ref. 11). The departure from isotropy of the particle veloc- 
ity distribution, caused by passage of the particles through successive apertures of the 
analyzer, was briefly considered by Boyd and Boylett (ref. 5). The effect of a magnetic 
field on the particle momentum distribution within an analyzer was discussed by Mason 
(ref. 4), by Simpson (ref. 8), and by Anderson, et al. (ref. 12). The integrated momen- 
tum distribution provided by straightforward application of the gridded analyzer is gener- 
ally deemed not as desirable a form of data presentation as the distribution function of the 
particles Itself. For this reason, means of electronically differentiating the integrated 
distribution have been developed by Boyd and Boylett (ref. 5), by Leder and Simpson 
(ref. 6), and by Krawec (ref. 13). 

These difficulties and the methods presented for dealing with data affected by them 
may be said to be well understood when the incident velocity distribution is isotropic in 
velocity space. There are many applications of the gridded electrostatic analyzer, how- 
ever, in which the relation between axial momentum and kinetic temperature is compli- 
cated by the fact that ions impinging on the analyzer are not isotropically distributed over 
a hemisphere in velocity space. 

This report discusses the results to be expected when the particles reaching the 
probe are "isotropic” within a cone in velocity space, as shown in figure 2. Such a ve- 
locity distribution may arise among the particles that escape through the ends of a "mag- 
netic mirror" containment geometry. Measurements made in such an experimental con- 
figuration have been reported by Eubank (ref. 1), Roth (ref. 2), and Zubov, et al. 

(ref. 14). The anal 3 dical form of the integrated momentum distribution does not seem to 
be available in the literature, except for a special case quoted without derivation by 
Eubank (ref. 1, p. 100). Analsrtical expressions are derived in the following analysis 
which relate the integrated spectrum of axial momentum observed by the collector to the 
effective escape cone angle, the floating potential, and the kinetic temperature of the par- 
ticles as they are incident on the analyzer. Because these curves are not simple expo- 


3 


Figure 2. - Schematic drawing of velocity space with 
particles confined to escape cone of half-angle 0^. 


nentials, a computer program (SPECTR) was written to obtain a best fit of the analytical 
expressions to the experimentally determined integrated spectrum. This program, writ- 
ten by Loretta R. Ellis, is presented in appendix B. 


ANALYTICAL CONSIDERATIONS 

Integrated Steady-State Spectra from Arbitrary Cone in Velocity Space 

In this section, it is assumed that the particles at the analyzer obey the Maxwellian 
distribution of velocities along a radius in velocity space and are "isotropic” within a 
cone in velocity space of half-angle 0^. It is also assumed that the velocity distribution 
is axisymmetric about cp, the azimuthal angle in velocity space. This is equivalent to 
the assumption that the particle motion in the magnetic field is "adiabatic" in the sense 
of preserving the magnetic moment of the particle (refs. 15 and 16), so that at a particu- 
lar location of the analyzer with respect to the plasma, the "escape cone" angle is 
not a function of particle energy. In order to make possible an analytical study of the in- 
tegrated momentum distribution, it is assumed that the magnetic field lines, electric 
field, and analyzer axis are parallel. The source of this plasma is assumed to float at 
some potential Vj with respect to the grounded zero reference of potential at which the 
analyzer is operated, so that an electric field (parallel to the analyzer axis and the mag- 
netic field) may exist between the plasma source and the analyzer. This floating potential 
is not that encountered in Langmuir probe theory. 

A particular physical situation represented by these assumptions (which represents 


Nonadiabatic 



CQ 



Figure 3. - Schematic drawing showing variation of escape cone angle 
of particles escaping a magnetic mirror confinement geometry. 


the experimental configuration actually used in testing the computer program) is shown in 
figure 3. The particles move along the magnetic field line and fall through a potential 
Vjj (the case of their acquiring energy is illustrated), so that all particles have veloci- 
ties along the z-axis greater than some minimum value. Because the particles behaved 
adiabatically in moving to the weaker magnetic field at position B, the escape cone angle 
is smaller than n/2. As the particles move to position C, the potential Vj2 ^ 
through which they fall is greater, and the minimum z-momentum increases further. The 
escape cone angle becomes still smaller in the much weaker magnetic field at C. It may 
be so weak that nonadiabatic effects occur (refs. 15 and 16), and the escape cone opens 
out above its adiabatic value at high particle energies. 

The plasma at the analyzer, according to the stated assiunptions, will have a distri- 
bution function in velocity space given by (ref. 17) 


dn(v, e>, (p) 




0 < 0 < 0 

o 


( 1 ) 
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where the symbols are defined in appendix A. The spatial number density at the analyzer 
will be n^ only if the distribution is completely isotropic (0^ = n), otherwise the density 
will be nQ(l - cos 0 q)/ 2. The most probable speed is given in terms of the kinetic 
temperature kT = eV^ by 


V 


o 



( 2 ) 


With the aid of figure 4, the speed v may be written in terms of the velocity v„ along 

z 

the analyzer axis, 


= V cos 0 

(3a) 

v_ sec 0 
z 

(3b) 


E equation (3b) is used to change the variable of integration from v to v^. in equa- 
tion (1), and an integration over the azimuthal angle 0 < ^ < 2ir is performed, one ob- 
tains 


dn(v2,, 0) 





sec 0 sin 0 dv„ 


d0 


( 4 ) 



Rgure4. - Geometric relation of velocities. 
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The differential current density in amperes per square meter of particles impinging 
on the analyzer is (ref. 17, pp. 206-209) 

0) = ev cos e dn(v, 0) = ev^, dn(v 2 ,, 9) (5) 


Substituting equation (4) into equation (5) results in 


2n e V, 

di(v„,0) = — - — exp|- 
^ .1/2 „3 


v^sec^O' 


I sec 0 tan 9 dv„ d0 


( 6 ) 


Integrating this velocity distribution over the cone in velocity space which contains the 
particles 0 ^ 0 0^ yields 



The current flowing to the collector plate may be determined by integrating equa- 
tion (7) over all v large enough to overcome the potential to which grid 2 is biased. 
z 

Since the particles may fall through the plasma floating potential Vj, the minimum ve- 
locity v^^.„ is 


'’^zmin 



( 8 ) 


where the plus sign is used if is positive, and the negative if Vj^ is negative. In 
practice, particles with negative v^, will never reach the analyzer, and the lower limit 
of integration is zero in this case. The analysis of data can be greatly simplified if 
Vj = 0. This can be accomplished in an experimental situation by biasing the analyzer 
ground to equal the floating potential Vj of the plasma source. 

The presence of a floating potential on the plasma source will result in translat- 
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ing the abscissa of the integrated energy distribution curve by an amount equal to Vj. 
When the analyzer is set at some potential Vj, the total transmitted current of posi- 
tive ions to the collecting electrode is given by integrating equation (7) over the range 

v„ < v„ < 00 , where 
zo z ’ 


v.o-V^<V-V.) 

and therefore the total ion current density in amperes per square meter is 


( 9 ) 



Performing the integration and substituting equations (2) and (9), one obtains, for the col- 
lector current as a function of retarding potential, 

r 


o o , 

] = — -T7 <(exp 


2jt 


1/2 


(V - Vj) 


- cos 0^ exp 


(V - Vj)sec^0^ 


( 11 ) 


V. 


J 


in amperes per square meter. The maximum possible current available to the analyzer 
occurs when V = Vj, because all particles emitted by the source are collected. This 
current is then given by 




( 12 ) 


If equation (12) is written in terms of the average speed v of the particles and if the par- 
ticles are assumed to be isotropically distributed over a hemisphere in velocity space, 
the result is 


h = 


n„ev 

o 

4 



( 13 ) 


This result for wall bombardment is familiar from kinetic theory (ref. 17, p. 208). 
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K equation (12) is used to normalize equation (11), one obtains 


3 = ij 


V< V, 


(14) 


-L = exp! 


(V- Vj) 





- 


1 - cos 6^ exp 

(V- 

Vj)tan20^ 


Vo 



sin^^o 


1 






VS V. 


(15) 


Equation (15) is the desired expression for the retarding potential curve. It is a 
function of the plasma floating potential the kinetic temperature of the plasma V^, 
and the effective escape cone angle in velocity space 6^. When Vj < 0, the grounded 
aperture and grid 3 will admit only the "tail" of the velocity distribution, since the low 
v_ particles will not have sufficient energy to overcome the potential barrier Vi and 
reach the analyzer. 

Consider some limiting cases of equation (15). When Vj « V^, or if the analyzer 
ground is biased to the plasma potential so that Vj = 0, equation (15) becomes 



VS 0 


(16) 


If, in addition, the cone angle 9^ is assumed to be determined by adiabatic motion of the 
particles from a magnetic mirror maximum to the magnetic field existing at the 

IllciX 

collector B^, the escape cone angle may be written as (ref. 18) 


Sln2e„=3- 

O T 5 


max 


(17) 


From equation (17) it follows that 


2q _ ®max " ®o 

cos 9^ = — 


B 


max 


(18) 
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and 


tan^e^ = 
o 


B “ B 
max o 


( 19 ) 


Substituting equations (18) and (19) into equation (16) results in the following equation for 
the normalized retarding potential curve: 



Equation (20) is identical with an expression for the collector current quoted without 
derivation by Eubank (ref. 1). From the preceding derivation, it is clear that equa- 
tion (20) holds for the particular case in which the analyzer ground is biased to the float- 
ing potential of the plasma, in which the particles move adiabatically along the magnetic 
field lines, and in which the particle distribution is uniform in angle and Maxwellian in 
velocity within the escape cone angle 0^ at the location of the probe. 

In the limit in which the particles are isotropically distributed over a hemisphere in 
velocity space, 0^ = tt/ 2, and equation (15) becomes simply 






( 21 ) 


As the escape cone angle becomes very small (as would happen if particles behaved adi- 
abatically in moving from a very strong magnetic field to a very weak field at the ana- 
lyzer), equation (15) becomes by L’Kopital’s rule (ref. 19, p. 591) 




h 



0, V> Vj 


( 22 ) 


A comparison of equations (21) and (22) reveals that narrowing tne escape cone angle en- 
hances the relative fraction of the particles above any given energy. 
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Integrated Pulsed Spectra from Arbitrary Cone in Velocity Space 


The preceding derivation could be described as assuming a steady state at the ana- 
lyzer in which the distribution function of the particles is zero outside a cone in velocity 
space and a Maxwellian function of velocity within it. In the subsequent analysis, a sec- 
ond model is examined in which the arrival rates, averaged over time, are proportional 
to the number of particles within each velocity increment dv. This situation could come 
about from the operation of a pulsed plasma source, such as laser- generated plasmas, 
theta pinches, and the experiment which is described in conjunction with this report. 

The limits of the cone in velocity space are determined by the geometry of the magnetic 
field in such applications. The assumption of a Maxwellian distribution, ’’isotropic” 
within the velocity space cone, at the site of the analyzer, may not be rigorously true in 
all experimental situations. However, this assumption is retained along with the other 
assumptions of the preceding case in order to simplify the analysis. The velocity- space 
distribution function of particles at the analyzer is obtained from equation (1) and (3b) 


dnCVg,, 9, (p) = 





sec^e sin 0 dv„ 


d0 d(p 


0 < 0 < 0^ (23) 


where the total number of particles reaching the analyzer per pulse is [N^(l - cos 0^)]/2. 
One cannot relate directly to a spatial number density at the analyzer because the 
particles with different velocities will be spread in space. If there are W bursts of par- 
ticles per second, the collected current in amperes will be 


di(v2, " 


2eN W v; 
o z 



9 9 

v^sec^e^x 9 

° 'sec'^0 tan 0 dv^, d0 


0 < 0^ (24) 


where the integration over the axis of symmetry (p has been performed. 

The current collected when the analyzer is set to repel all particles within the escape 

cone 0„ and with a z-momentum less than that corresponding to v„. is given by 
o zo 



11 



Performing the integration over the angle 0, one obtains 


1 = 


eN W 
o 



- exp 


v^sec^e 
z '^o 


zo 


dv„ 


(26) 


Performing the integration over v and making use of equations (2) and (9) result in 


eN„W 
i = — 


erfcl 


(V- V. 


1/2 




- cos 0^ erfc 


(V- Vj)sec\ 


1/2 


V> V. 


(27) 


The maximum possible current to the analyzer occurs when V = Vj and is given by 


eN W 

il=-^(l -COS0^) 


(28) 


Using equation (28) to normalize equation (27) gives the shape of the retarding potential 


curve: 


1 = ii 


V< V, 


(29) 


and 


1/2 


V> Vj (30) 

The identity erfc x s 1 - erf x has been used in equation (30) . 

Some limiting cases of equation (30) will be examined. When the velocity distribution 
is isotropic over a hemisphere in velocity space, and 0^ = tt/2, equation (30) becomes 


erfl 


(V- V. 


1/2 


- cos 0^ erf 


-L=l 


(V - Vj) 2 
^ sec^0. 


V, 


1 - cos 0, 
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— = 1 - erf 


(Y- V, 


. 1/2 


Vi, 0^=^ 


(31) 


On the contrary, when the particles are confined within a very narrow cone in velocity 
space, equation (30) becomes, by L’Hopital’s rule. 


r 

‘1 


x2e-=^ dx 


va Vj, 


(32) 


where 




(33) 


Other Integrated Energy Spectra 

In a particular physical situation, arguments may be available which would lead to 
prediction of particular types of anisotropic and non-Maxwellian distributions. Calcula- 
tions similar to those presented for the two preceding models should then be performed. 
The pertinent details of the plasma processes usually are not well enough understood to 
justify the choice of an elaborate model. The Maxwelllanized plasmas discussed pre- 
viously may be considered the limiting case of a perfectly random distribution. It is rec- 
ognized that, in the first case, the assumption of a Maxwellian distribution at the ana- 
lyzer implies a slightly non-Maxwellian distribution at other axial locations, if adiabatic 
motion is assumed. The two cases considered are in a way extremes - they give the 
maximum and minimum arrival rates (and hence V^) consistent with locally Maxwellian 
distributions. 

It is often of interest to compare these two cases with the opposite limit of a per- 
fectly ordered distribution. Two such ordered distributions were included in the com- 
puter program. The first is the delta- function energy distribution shown in figure 5(a). 
The integrated energy distribution shown in figure 5(b) would be a square wave of the 
form 


1= U 


V< 


(34) 
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(a) Delta-function energy distribution. 
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(b) Integrated delta-function distribution. 

Figure 5. - Schematic drawing of delta- 
function energy spectra. 
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(a) Square-wave distribution function. 



(b) Integrated square-wave distribution. 

Figure 6. - Schematic drawing of square- 
wave energy spectra. 


and 


i = 0 


V> V. 


( 35 ) 


This distribution would be appropriate to a perfectly monoenergetic beam of energy 
(in electron volts), whose velocity is directed along the analyzer axis. 

A second simple distribution function is the "square wave" distribution, in which all 
the particles reaching the analyzer are evenly distributed in energy between ^ V V 2 , 
as illustrated in figure 6(a). The particles have only a negligible spread in angle in ve- 
locity space about the analyzer axis. Such a particle distribution will yield the integrated 
distribution shown in figure 6(b), whose anal 5 rtic form is 


V< Vj 


( 36 ) 


i = 



Vj ^ V< 


Vo 


( 37 ) 
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i = 0 


(38) 


V>V2 

An average ’’temperature” for this distribution can be defined by 


(39) 


COMPUTER PROGRAM FOR DATA ANALYSIS 


Equations (15) and (30) for the integrated energy spectra are smooth, monotonically 
decreasing functions of the retarding potential V. Unfortunately, experimental data 
obeying these analytical relations cannot be reduced in any simple way, since the curves 
are not, in general, simple exponentials. In the special case of an isotropic steady- state 
distribution over the entire hemisphere in velocity space, equation (21) shows that the in- 
tegrated distribution is exponential. The kinetic temperature of the plasma in question 
may then be determined by plotting the data on a semilogarithmic graph and measuring 
the slope of the resulting straight line. When 9^ * ir/2, however, the integrated energy 
spectra can depart significantly from an exponential curve, and an attempt to force an ex- 
ponential fit may result in erroneous values of the kinetic temperature V^, and/or the 
floating potential Vj. 

A computer program was written which provides an iterated best fit of the four inte- 
grated energy spectra discussed previously to a tabulated set of experimental data. This 
computer program, its operation, and its application to data analysis are discussed in ap- 
pendix B. The integrated distribution functions actually written into the computer pro- 
gram are as follows: In the steady- state case, from equation (15), 


1 = ^0 -^II 



= Io-^Il 

v< 




(V - vjtan^e 

1 - cos 0Q exp 

^ 


o 

L 

J 


For the pulsed case, from equation (30), 


(40) 


VS Vj (41) 


I = Io + Ii V<Vi 


(42) 
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r 


erfl 


rv- V. 


.1/2 


- cos 0Q erfl 


I = Io + Ii/ 


1 - 


(V - Vj)sec^0^ 


1/2 


1 - cos Q, 




For the delta- function distribution, from equations (34) and (35), 


I = Io-Mi V<V^ 


V ^ Vj (43) 


(44) 

(45) 


For the "square-wave" distribution, from equations (36) to (38), 

I = + Ij V < Vj (46) 

/V2 - V\ 

1 = In + ll V. < V < V„ (47) 

° \ V 2 -V 

1 = 1^ V > ^2 

In the preceding equations, the parameter represents the zero-offset of the curve on 
the y-axis and can be made less than a few percent of Ij by zeroing the X-Y recorder 
before each run. 

The computer program takes as input the retarding potential curves in the form of a 
table of values of V and its corresponding collector current I. Through an iterative 
process, the program converges to best-fitting values of the five parameters I^, Ij, V^, 
Vj, and 0^, appropriate to the distribution function in question. The computer output 
illustrated in table I consists of the raw data for V and I in the first two columns, and 
the best-fitting curve points for the four distributions given in equations (40) to (48) in the 
next four columns. Each row of the output table represents one value of this retarding 
potential at which a point was tabulated from the original X-Y recording. At the bottom 
of each column the best fitting values of I^, 1^, V^, V^, and 0^ for that distribution are 
printed out. The mean square error of the best-fitting curve is also printed out. 

The parameter represents the total flux of particles reaching the analyzer, when 
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TABLE I. - EXAMPLE OF OUTPUT FROM COMPUTER PROGRAM 


[V is in volts and I is 


RUN NUMBER AS- 3 



V 

I 


STEADY 

1 

0. 

0.14100E 

03 

0.14121E 03 

2 

50.000 

0.13400E 

03 

0.12784E 03 

3 

100.000 

0.11780E 

03 

0.11570E 03 

4 

150.000 

0.10200E 

03 

0.10471E 03 

5 

200.000 

0.92300E 

02 

0.94746E 02 

6 

250.000 

0.82500E 

02 

0.85721E 02 

7 

300.000 

0.75000E 

02 

0.77547E 02 

8 

350.000 

0.69000E 

02 

0.70142E 02 

9 

400.000 

0.63000E 

02 

0.63435E 02 

10 

450.000 

0.56600E 

02 

0.57359E 02 

11 

500.000 

0.51750E 

02 

0.51856E 02 

12 

550.000 

0.46500E 

02 

0.46871E 02 

13 

600.000 

0.44250E 

02 

0.42356E 02 

14 

650.000 

0.37500E 

02 

0.38265E 02 

15 

700.000 

0.33000E 

02 

0.34560E 02 

16 

750.000 

0.30000E 

02 

0.31205E 02 

17 

800.000 

0.27000E 

02 

0.28165E 02 

18 

850.000 

0.25000E 

02 

0.25411E 02 

19 

900.000 

0.22500E 

02 

0.22917E 02 

20 

950.000 

0.21000E 

02 

0.20658E 02 

21 

1000.000 

0.19000E 

02 

0.18611E 02 

22 

1050.000 

0.17000E 

02 

0.16757E 02 

23 

1100.000 

0.16000E 

02 

0.15078E 02 

24 

1150.000 

0.14500E 

02 

0.13557E 02 

25 

1200.000 

0.13500E 

02 

0.12180E 02 

26 

1250.000 

0.12500E 

02 

0.10932E 02 

27 

1300.000 

O-llOOOE 

02 

0.98012E 01 

28 

1350.000 

0.10500E 

02 

0.87773E 01 

29 

1400.000 

O.IOOOOE 

02 

0.78498E 01 

30 

1450.000 

0.90000E 

01 

0.70097E 01 

31 

1500.000 

0.85000E 

01 

0.62487E 01 

32 

1550.000 

0.75000E 

01 

0.55593E 01 

33 

1600.000 

0.70000E 

01 

0.49349E 01 

34 

1650.000 

0.60000E 

01 

0.43693E 01 

35 

1700.000 

0-55000E 

01 

0.38570E 01 

36 

1750.000 

0.50000E 

01 

0.33930E 01 

37 

1800.000 

0.45000E 

01 

0.29726E 01 

38 

1850.000 

0.40000E 

01 

0.25918E 01 

39 

1900.000 

0.35000E 

01 

0.22469E 01 

40 

1950.000 

0.30000E 

01 

0.19345E 01 

41 

2000.000 

0.30000E 

01 

0.16515E 01 

42 

2050.000 

0.20000E 

01 

0.13952E 01 

43 

2100.000 

O.IOOOOE 

01 

0.11630E 01 

44 

2150.000 

O.IOOOOE 

01 

0.95268E 00 

45 

2200.000 

0.50000E 

00 

0.76217E 00 

46 

2250.000 

0. 


0.58960E 00 

47 

2300.000 

0. 


0.43329E 00 

48 

2350.000 

-0.50000E 

00 

0.29170E 00 

49 

2400.000 

-O.IOOOOE 

01 

0.16344E 00 

50 

2450.000 

-O.IOOOOE 

01 

0.47266E-01 

51 

2500.000 

-0.15000E 

01 

-0.57967E-01 



INITIAL 


STEADY 


10 = 

-0.80000E 

01 

-0.10700E 

01 

11 = 

0.14500E 

03 

0.22431E 

03 

VO = 

0.35000E 

03 

0.50544E 

03 

VI = 

0.30000E 

02 

-0.27732E 

03 

THETA = 

0.40000E 

02 

0.72591E 

02 

ERROR = 



0.12815E 

01 


amperes xlO . ] 


SQ-WAVE DELTA PULSED 

0,12154E 03 0.89541E 02 0.14417E 03 
0.11544E 03 0.89541E 02 0.12897E 03 
0.10934E 03 0.89541E 02 0,11559E 03 
0.10324E 03 0.89541E 02 0.10381E 03 
0.97137E 02 0.89541E 02 0.93402E 02 
0.91037E 02 0.89541E 02 0.84189E 02 
0.84937E 02 0.89541E 02 0.76004E 02 
0.78837E 02 0.89541E 02 0.68710E 02 
0.72737E 02 0.89541E 02 0.62188E 02 
0.66637E 02 0.89541E 02 0.56342E 02 
0.60537E 02 0.99946E 01 0.51088E 02 
0.54437E 02 0.99946E 01 0.46355E 02 
0.48337E 02 0.99946E 01 0.42083E 02 
0.42237E 02 0.99946E 01 0.38220E 02 
0.36137E 02 0.99946E 01 0.34721E 02 
0.30037E 02 0.99946E 01 0.31547E 02 
0.23937E 02 0.99946E 01 0.28664E 02 
0.17837E 02 0.99946E 01 0.26042E 02 
0.11737E 02 0.99946E 01 0.23656E 02 
0.56365E 01 0.99946E 01 0.21480E 02 
0.46745E 01 0.99946E 01 0.19496E 02 
0.46745E 01 0.99946E 01 0.17685E 02 
0.46745E 01 0.99946E 01 0.16030E 02 
0.46745E 01 0.99946E 01 0.14517E 02 
0.46745E 01 0.99946E 01 0.13132E 02 
0.46745E 01 0.99946E 01 0.11865E 02 
0.46745E 01 0.99946E 01 0.10703E 02 
0.46745E 01 0-99946E 01 0.96392E 01 
0.46745E 01 0.99946E 01 0.86632E 01 
0.46745E 01 0.99946E 01 0.77677E 01 
0.46745E 01 0.99946E 01 0.69458E 01 
0.46745E 01 0.99946E 01 0.61911E 01 
0.46745E 01 0.99946E 01 0.54977E 01 
0.46745E 01 0.99946E 01 0.48605E 01 
0.46745E 01 0.99946E 01 0.42747E 01 
0.46745E 01 0.99946E 01 0.37360E 01 
0.46745E 01 0.99946E 01 0.32404E 01 
0.46745E 01 0.99946E 01 0.27843E 01 
0.46745E 01 0.99946E 01 0.23645E 01 
0.46745E 01 0.99946E 01 0.19780E 01 
0.46745E 01 0.99946E 01 0.16219E 01 
0.46745E 01 0.99946E 01 0.12940E 01 
0.46745E 01 0.99946E 01 0.99174E 00 
0.46745E 01 0.99946E 01 0.71319E 00 
0.46745E 01 0.99946E 01 0.45640E 00 
0.46745E 01 0.99946E 01 0-21961E 00 
0.46745E 01 0.99946E 01 0.12210E-02 
0.46745E 01 0.99946E 01 -0.20023E 00 
0.46745E 01 0.99946E 01 -0.38610E 00 
0.46745E 01 0.99946E 01 -0.55762E 00 
0.46745E 01 0.99946E 01 -0.71593E 00 


SQ-WAVE DELTA PULSED 

0.46745E 01 0.99946E 01 -0.26540E 01 

0.12403E 03 0.79546E 02 0.23439E 03 

0.50833E 03 0.50000E 03 0.70571E 03 

-0.58784E 02 0. -0.23686E 03 

0. 0. 0.69349E 02 

0.54365E 01 0.12354E 02 0.10056E 01 
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the retarding potential is V = V^. In the case of a steady-state flux, can be written 
by using equation (12), 


n^ev A « 

I. = -2_o_ G sin^0 

1 I/O O 


27T 


1/2 


and by using equation (28) for the pulsed case. 


eN W 

h = ^o) 


(49) 


(50) 


The area of the analyzer opening is A, and the factor G in equation (49) accounts for 
geometrical effects in the analyzer, such as the transmissivity of the grids. The velocity 
v^ can be calculated from the kinetic temperature of the plasma, and equation (49) 
then gives Ij as proportional to the number density n^ of the plasma source. The pa- 
rameter Vj is the effective floating potential of the plasma source at the particular field 
line on which the analyzer is located, and 9 ^ is the effective escape cone angle of the 
particles as they enter the analyzer. It should always be kept in mind that the best- 
fitting values of 1^, Ij, V^, Vj, and 9 ^ are effective values only. These values can be 
identified with the actual properties of the plasma only when the physical model is appli- 
cable and when the systematic errors that might affect them are negligible. 

The computer program requires starting values for I^, 1^, V^, and Vj (a universal 
starting value of 0^ = 40° is written into the program). It is not necessary to make very 
precise estimates of these quantities, since the computer program will converge to the 
best-fitting values even if the initial estimates are off by a factor of 2 or more. It is not 
difficult with a little practice to examine an X-Y plot of the retarding I against V curve, 
and estimate 1^ and 1^ by inspection to within 5 percent, and and Vj to within a 
factor of 2, when Vj is positive. The computer time required for an iterated fit to a set 
of experimental data was about 6 seconds on an IBM 709411/7044 direct couple computer 
system. 

Experience showed that at least 12 points should be read along the I against V re- 
tarding potential curves in order for the computer program to obtain a good fit. The tab- 
ulated data points which form the input to the program should be taken at evenly spaced 
intervals of V along the entire retarding potential curve. If the data points were bunched 
in a particular region of the retarding- potential curve, the computer program would 
weigh that region proportionately more heavily than regions with only a few data points. 

The mean square error is based on the difference of the entire best-fitting curve and 
the raw data normalized to the value of I when V = 0. A mean square error of 1 per- 
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cent, for example, means that the average data point deviates from the best-fitting curve 
by an amount that is 1 percent of I when the retarding potential is zero. Plotting a large 
number of best-fitting curves on the corresponding I - V data plots showed that a mean 
square error of less than 1 percent was a very good fit, 1 to 1. 5 percent was a good fit, 

1. 5 to 2. 0 percent was only a fair fit, and greater than 2 percent was a relatively poor fit. 
As the mean square error becomes progressively greater than 1. 5 percent, it becomes 
progressively more probable that either the velocity distribution is non-Maxwellian, or 
that some systematic error has affected the retarding-potential curve. 


EXPERIMENTAL APPARATUS AND PROCEDURE 
Description of Apparatus 

The superconducting magnet facility in which the experimental data were taken has 
been described previously (ref. 20). The properties of the plasma and the means used to 
generate it have been described elsewhere (refs. 2 and 21). The discharge characteris- 
tics depend on anode voltage, vacuum tank pressure, etc. Of relevance to the present 
discussion is the observation of strong turbulence in the plasma, which apparently pro- 
vides the means of randomizing the velocity distribution, increasing the kinetic tempera- 
ture through turbulent heating, and scattering particles from this randomized distribution 
into the escape cone (ref. 21). 

A scale drawing of the analyzer used in the experiments to be described subsequently 
is shown in figure 1. The circuitry and location of the analyzer are shown in figure 7. 
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Rgure 7. - Schematic drawing of circuitry for electrostatic analyzer. 
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Illllllllll 


The analyzer is mounted on a movable support which enables it to be displaced in a direc- 
tion perpendicular to the axis of the magnetic field in a plane containing the magnetic field 
axis. An X-Y recorder was used to record the retarding potential curves in the present 
experiment. A fast- sweep oscilloscope display would be desirable in an experiment of 
short duration. 


Experimental Procedure 

The front surface of the analyzer was placed 20 centimeters away from the point of 
maximum magnetic field. The axis of the probe was parallel to the magnetic field axis 
but displaced radially 10 centimeters from it. This placement was deemed desirable be- 
cause the ion current impinging on the analyzer near the magnetic field axis was too noisy 
to produce smooth retarding potential curves of high quality. Unfortunately this place- 
ment resulted in the introduction of a systematic error in the effective escape cone angle 
e^, since the analyzer axis was at an angle of about 30° to the local magnetic field lines. 

The field strength at the magnetic field maximum was 1. 0 tesla, about 0. 15 tesla at 
the front surface of the analyzer and about 0. 10 tesla at the collector. The magnetic 
fields at the location of the energy probe were small enough to cause nonadiabatic effects 
among ions of the kinetic temperatures observed (refs. 15 and 16), and one would there- 
fore not necessarily expect ions of high energy to be confined to an escape cone whose 
angle is given by equation (17). 

The plasma was observed to float at a potential that is typically from 
-300 :< Vj ^ 200 volts. The electron energies within the plasma appeared to be about 
equal to the floating potential. 

The collector currents were in all cases below 10"^ ampere, and the corresponding 
number density of particles in the beam was below 10° particles per cubic centimeter for 
the results presented herein. The minimum electron Debye length was larger than the 
analyzer diameter, and the local value of the ion gyroradius, corresponding to the ion 
kinetic temperature, was much larger than the diameter of the grid holes. 

During operation, the 2. 54-centimeter-diameter aperture and grid 3 were grounded, 
as shown on figure 7. Grid 3 was used in preliminary measurements to verify that the 
electron current reaching the analyzer was much smaller than the ion current. Grid 2 
was connected to a well-filtered 5-kilovolt d-c power supply. This power supply could be 
varied from zero to full scale by a motor drive and variable transformer in the power 
supply primary circuit. This assured that the retarding potential appearing on grid 2 
would vary smoothly and continuously in time. Grid 1 was biased from -200 to -600 volts 
dc. Its function was the repelling of secondary electrons that might be emitted by ions 
that struck the collector, and also the repelling of any electrons below this energy that 
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entered with the ion beam. The collector was a flat brass disk that was connected to an 
electronic microvolt-ammeter. The output of this instrument was fed into the Y-axis of 
an X-Y recorder. The X-axis of the recorder was driven by a voltage proportional to 
the retarding potential applied on grid 2. Approximately 1 minute was required to sweep 
out a complete retarding- potential curve. 

The X-Y graphs of the retarding potential curves, illustrated in figure 8, were con- 
verted to tabular form by picking and tabulating approximately 50 values of the collector 
current for evenly spaced values of the retarding potential. An estimate was then made 
of the values of 1^, Ij, V^, and Vj to start the iterative process. These data were 
then placed on punched cards and used as input to the computer program. 


EXPERIMENTAL RESULTS 
Illustrative Experimental Data 

Four series of experimental runs, whose characteristics are summarized in table II, 
were taken with the energy analyzer. A run consisted of taking a single retarding poten- 
tial curve; a series consisted of many runs taken under the same set of operating condi- 
tions. These series of runs illustrate the application of the computer program to an ac- 
tual set of experimental data and make it possible to investigate the magnitude of the sys- 
tematic and random errors encountered in practice. 

The second column of table II lists the number of runs taken in each series. Each 
series listed in table n was made on the same day, in sequence, and under conditions 
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TABLE n. - OPERATING PARAMETERS FOR RUN SERIES INVESTIGATED 


Run 

series 

Number 
of runs 
in series 

Gas used 

Tank pi 
torr 

fessure 

N/m^ 

Bias 
voltage 
on grid 1, 
V 

Anode 

current, 

A 

Elapsed 
time for 
series, 
min 

x-axis scale, 
retarding 
potential on 
grid 2, 
V/in. 

y-axis scale, 
collector 
current. 
A/in. 

AS 

10 

Helium 

9 X10‘® 

1. 2xl0‘® 

-200 

0. 005 

76 

200 

axio"*^ 

AT 

10 

Helium 

4.9x10"® 

6. 5x10'® 

-500 

. 061 

36 

200 

10-® 

AU 

13 

Neon 

1. 9xl0‘® 

2. 5x10"'^ 

-500 

.010 

43 

200 

lO""^ 

AV 

13 

Neon 

4. 0xl0‘® 

5. 3xl0“* 

-500 

.030 

40 

200 

axio""^ 

CK-CL 

26 

Deuterium 

9. 0x10"® 

1. 2x10"® 

-600 

. 004 

— 

500 

10"'^ 

CM-CN 

17 

Deuterium 

1. 2x10" ® 

1. 6x10" ® 

-600 

1 .0025 

— 

200 

10""^ 


(vacuum tank pressure, plasma number density, etc. , ) that were as nearly identical as 
the nature and limitations of the apparatus would permit. The duration over which a 
series was taken is listed in the seventh column of table n. The quantity that was ob- 
served to vary most during a series of runs was the vacuum tank pressure. The bias 
voltage used on grid 1 to repel secondary electrons emitted from the collector is listed 
in column 5. The scales of the ordinate and abscissa of the retarding potential curves 
from the X-Y recorder are shown in the last two columns of table n. A typical X-Y plot 
of the retarding potential from the AS series of runs is shown in figure 8. 


Computer Generated Best-Fitting Integrated Spectra 

The experimental retarding potential curves described for the first four series of 
runs listed in table n (altogether 46 curves) were tabulated and used as input to the com- 
puter program, which then obtained best-fitting theoretical curves. The computer output 
for the experimental run AS- 3, shown in figure 8, is given in table I. Data were read from 
the curve in figure 8 for 51 values of the retarding potential V. These raw data points of 
the retarding potential V and collector current I then appear in the first two columns of 
table I. From 12 to 20 data points are sufficient to define the shape of the curve for use 
by the computer program. Between 50 and 60 raw data points were used for each run 
listed in table II, since these runs were being used for illustrative purposes and also be- 
cause it was desired to eliminate an inadequate number of points as a possible source of 
error. 

The last four columns in table I are the analsrtical curves that best fit the experi- 
mental data for the four possible integrated energy spectra discussed previously. At the 
bottom of table I are shown the best-fitting values of 1^, Ij, V^, Vj, and 6 ^ for the four 
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TABLE m. - ANALYSIS OF AS RUN SERIES 


(a) Parameters for best-fitting steady curve 


AS series 

y-axis zero 

Total particle 

Kinetic 

Floating 

Effective 

Mean square 

Elapsed time, 

Tank pressure 

run number 

offset, 

flux reaching 

temper- 

potential, 

escape cone 

error. 

min 

. 



2 


Io» 

analyzer, 

ature, 

Vp 

angle in 

percent 


torr 

iN/m 


A 

h> 


V 

velocity 









A 

eV/ion 


space, 












®o» 

deg 







3 

-1. 07x10"® 

224x10"® 

505 

-277 

73 

1. 28 

9 

9. IxlO"® 

1.21x10"® 

6 

-1. 53 

218 

501 

-303 

69 

1. 49 

22 





9 

36 

268 

508 

-355 

85 

1. 36 

33 





12 

-1. 60 

251 

510 

-325 

83 

1. 46 

43 





15 

-1. 14 

215 

513 

-335 

67 

1. 36 

51 

9.0x10"® 

1. 195X10"® 

18 

-.48 

235 

506 

-345 

73 

1. 53 

59 

9.0 


1. 195 

20 

77 

244 

526 

-362 

78 

1. 37 

64 

8. 85 

1. 175 

22 

67 

286 

519 

-446 

76 

1. 49 

68 

8. 80 


1. 17 


24 

-1. 07 

229 

515 

-332 

72 

1.73 

72 

9.0 


1. 195 

26 

.65 

231 

493 

-326 

72 

1. 29 

76 

9. 05 

1. 20 


Average 

-0. 80x10"® 

240x10"® 

510 

-341 

75 

1. 44 

— 

9. OlxlO"® 

1. 196x10"® 

Mean error 

±0. 48x10"® 

±18x10"® 

±7 

±19 

±5 

±0. 10 

- 

0. 08x10"® 

0. 016x10"® 


(b) Parameters for best-fitting pulsed curve 


AS series 

y-axis zero offset, 

Total particle 

Kinetic 

Floating potential. 

Effective escape cone 

Mean square error, 

run number 

^o’ 

flux reaching 

temperature, 

Vp 

angle in velocity space, 

percent 


A 

analyzer, 

Vo. 

V 

00. 




Il« 

eV/ion 


deg 




A 





3 

-2. 65x10"® 

234X10"® 

706 

-237 

69 

1. 01 

6 

-3. 04 

339 

697 

-265 

85 

1. 04 

9 

-2. 21 

223 

736 

-198 

75 

. 85 

12 

-3. 09 

331 

708 

-275 

84 

1. 08 

15 

-2. 58 

297 

709 

-291 

78 

. 98 

18 

-1.79 

341 

693 

-297 

84 

1. 13 

20 

-2. 10 

357 

718 

-314 

87 

. 96 

22 

-2, 03 

366 

710 

-308 

88 

1. 11 

24 

-2, 49 

238 

714 

-266 

71 

1. 41 

26 

-. 80 

290 

681 

-279 

77 

. 92 

Average 

-2. 28X10"® 

302x10"® 

707 

-273 

80 

1, 05 

Mean error 

±0. 49x10"® 

±42x10"® 

±10 

±25 

±6 

±0. 12 


(c) Parameters for best- fitting square-wave curve (d) Parameters for best- fitting delta- function curve 


AS series 

y-axis zero 

Total particle 

Kinetic 

Floating 

Mean square 


AS series 

y-axis zero 

Total particle 

Kinetic 

Mean square 

run number 

offset, 

flux reaching 

temper- 

potential, 

error, 


run number 

offset, 

flux reaching 

temper- 

error. 


lo- 

analyzer. 

atvLre, 

Vp 

percent 



lo- 

analyzer, 

ature, 

percent 


A 


Vo. 

V 




A 

Ip 

Vo- 




A 

eV/ion 






A 

eV/ion 


3 

4. 67x10"® 

124x10"® 

508 

-59 

5. 44 


3 

10. 0 XIO"® 

79. 5X10"® 

500 

12.4 

6 

3. 87 

125 

534 

-111 

5.61 


6 

8. 92 

76. 0 



12. 1 

9 

5. 17 

122 

532 

-106 

5. 53 


9 

10. 2 

74.9 



12. 2 

12 

3. 98 

124 

540 

-119 

5. 56 


12 

9. 41 

73. 5 



12.4 

15 

4. 42 

123 

549 

-134 

5. 56 


15 

9.08 

75. 2 



12.3 

18 

4. 82 

122 

550 

-142 

5. 58 


18 

9. 67 

72. 2 



12.1 

20 

4. 66 

122 

584 

-168 

5. 58 


20 

10.0 

72. 1 



12. 2 

22 

4. 99 

122 

554 

-144 

5. 68 


22 

9. 90 

72. 1 



12. 3 

24 

4.62 

123 

543 

-122 

5. 83 


24 

9.72 

74. 1 



12.6 

26 

6. 43 

120 

506 

-97 

5. 47 


26 

11.9 

74.9 

450 

12. 3 

Average 

4. 76x10"® 

123X10"® 

540 

-120 

5. 58 


Average 

9. 88X10"® 

74. 5X10"® 

495 

12. 3 

Mean error 

±0. 35X10"® 

±1. IXlO"® 

±16 

±22 

±0. 07 


Mean error 

±0. 43x10"® 

±1. 7x10"® 

±9 

±0.1 




TABLE IV. ~ ANALYSIS OF AT RUN SERIES 


(a) Parameters for best-fitting steady curve 


AT series 

y-axis zero 

Total particle 

Kinetic 

Floating 

Effective 

Mean square 

Elapsed time, 

Tank pressure | 

run number 

offset. 

flux reaching 

temper- 

potential, 

escape cone 

error, 

min 

. 


.. 

2 


Io> 

analyzer, 

ature. 

vi- 

angle in 

percent 


torr 

N/m 


A 

h> 

Vo» 

V 

velocity 









A 

eV/ion 


space, 












®0’ 












deg 







1 

-13. 9X10"® 

812x10"® 

391 

-46 

44 

0.75 

0 

4. 95x10"® 

6. 60X10"® 

3 

-12. 5 

784 

367 

-43 

41 

.77 

4 





5 

-10.4 

789 

365 

-35 

42 

.78 

8 





7 

-9.9 

782 

367 

-37 

42 

.76 

12 





9 

-10.7 

773 

358 

-55 

40 

.82 

16 





11 

-10.2 

760 

347 

-48 

38 

. 87 

20 

4. 90X10"® 

6.52x10"® 

13 

-7. 9 

754 

335 

-39 

38 

. 82 

24 

4.95 


6. 60 


15 

-15.0 

813 

401 

-63 

45 

.82 

27 

4. 90 


8. 52 


17 

-14.8 

794 

392 

-60 

43 

.98 

31 

4. 90 


6. 52 


19 

-17.6 

820 

422 

-65 

45 

.90 

34 

4. 90 


6. 52 


Average 

-12. 3X10"® 

788X10"® 

375 

-49 

42 

0. 83 

- 

4. 93x10"® 

6. 56X10"® 

Mean error 

±2. &<10"® 

±18x10"® 

±22 

±9 

±2 

±0. 06 

- 

±0. 0^10"® 

0. 026x10"® 


(b) Parameters for best-fitting pulsed curve 


AT series 

y-axis zero offset, 

Total particle 

Kinetic 

Floating potential, 

Effective escape cone 

Mean square error, 

run number 


flux reaching 

temperature, 

Vp 

angle in velocity space, 

percent 


A 

analyzer, 

Vo» 

V 

^o» 




Ip 

eV/ion 


deg 




A 





1 

-17. 8x10'® 

792X10"® 

464 

33 

41 

0. 53 

3 

-16.2 

765 

435 

36 

39 

. 51 

5 

-14.4 

772 

435 

40 

39 

.45 

7 

-14.0 

764 

437 

40 

39 

.44 

9 

-14.0 

749 

419 

29 

37 

. 54 

11 

-13.4 

737 

409 

35 

36 

. 56 

13 

-11.3 

734 

396 

39 

36 

.44 

15 

-18.2 

785 

468 

25 

41 

. 58 

17 

-18.2 

768 

458 

29 

39 

.75 

19 

-21. 0 

791 

493 

27 

41 

.67 

Average 

-15. 9X10"® 

p 

766x10 ^ 

441 

33 

39 

0. 55 

Mean error 

±2. 4x10"® 

±16x10"® 

±23 

±5 

±1 

±0. 08 


(c) Parameters for best-fitting square-wave curve (d) Parameters for best-fitting delta-fiinction curve 


AT series 

y-axis zero 

Total particle 

Kinetic 

Floating 

Mean square 

AT series 

y-axis zero 

Total particle 

Kinetic 

Mean square 

run number 

offset. 

flux reaching 

temper- 

potential, 

error, 

run number 

offset. 

flux reaching 

temper- 

error. 


V 

analyzer, 

ature. 

Vp 

percent 


Io» 

analyzer, 

ature. 

percent 


A 

Ip 

Vo. 

V 



A 

h’ 

Vo» 




A 

eV/ion 





A 

eV/ion 


1 

18. 0X10"® 

741x10"® 

509 

-8.0 

4. 28 

1 

68. 5x10"® 

501x10"® 

500 

14.2 

3 

17.6 

740 

502 

-31.2 

4. 28 

3 

61.9 

486 



14.1 

5 

19.0 

738 

493 

-16.7 

4.38 

5 

63. 4 

492 



14.1 

7 

17. 5 

740 

518 

-36.1 

4. 40 

7 

65. 1 

488 



14.2 

9 

16.3 

741 

513 

-62.4 

4. 36 

9 

58.1 

470 



14.0 

11 

15.8 

742 

516 

-71.7 

4. 35 

11 

56.7 

467 



14.0 

13 

15.4 

742 

512 

-73.1 

4. 42 

13 

61.7 

482 

450 

14.1 

15 

18.4 

739 

519 

-25.6 

4. 42 

15 

68.6 

492 

500 

14.3 

17 

17.3 

740 

529 

-43.2 

4. 40 

17 

67.2 

484 

500 

14. 5 

19 

18.2 

740 

540 

-22.8 

4.45 

19 

67.6 

483 

550 

14.4 

... 

Average 

17. 4x10"® 

740X10"® 

515 

-39.1 

4.37 

Average 

63. 9x10"® 

48^10"® 

500 

14.2 

Mean error 

±0. 9X10"® 

±1x10"® 

±9 

±18.8 

±0.05 

Mean error 

±3. 5x10"® 

±7x10"® 

±10 

±0.1 






TABLE V. - ANALYSIS OF AU RUN SERIES 


(a) Parameters for best- fitting steady curve 


AU series 

y-axls zero 

Total particle 

Kinetic 

Floating 

Effective 

Mean square 

Elapsed time, 

Tank pressure 

run number 

offset, 

Io» 

A 

flux reaching 
analyzer, 

h’ 

A 

temper- 

ature, 

Vo. 

eV/ion 

potential, 

Vl. 

V 

escape cone 
angle in 
velocity 
space, 

00- 

deg 

error, 

percent 

min 

torr 

N/m^ 

1 

-1. 31x10"® 

47. 6xl0“® 

683 

-326 

43 

1. 38 

0 

1.7x10"® 

2. 26x10"® 

3 

-1. 10 

58.7 

662 

-199 

60 

2. 37 

7 

2. 0 

2. 66 

5 

47 

50.6 

592 

-154 

49 

.79 

9 

1. 95 

2, 60 

7 

-.80 

55. 9 

624 

-204 

56 

.71 

13 

2. 1 

2. 80 

9 

-1.17 

51,6 

703 

-387 

50 

1. 15 

17 

1.6 

2, 13 

11 

-1. 58 

53. 5 

727 

-354 

51 

1. 52 

20 

1.8 

2. 40 

13 

-. 40 

44. 3 

530 

-70 

38 

1. 11 

23 

1. 95 

2. 60 

15 

84 

48.1 

620 

-237 

45 

.79 

27 

1.85 

2. 46 

17 

83 

50. 1 

653 

-271 

48 

.71 

29 

1. 90 

2. 53 

19 

72 

53.7 

632 

-167 

53 

. 87 

35 

2.6 

3. 46 

21 

-1. 20 

42. 5 

667 

-331 

35 

.77 

39 

1.45 

1. 93 

22 

-1.44 

43. 8 

696 

-394 

37 

1. 17 

41 

1. 50 

2. 00 

23 

-. 91 

46. 0 

659 

-317 

41 

. 99 

43 

1. 55 

2. 06 

Average 

-0. 98X10“® 

49. 7X10"® 

650 

-262 

47 

1.10 

- 

1. 84x10"® 

2. 45x10"^ 

Mean error 

±0. 29x10"® 

±4. 0X10“® 

±39 

±84 

±6 

±0. 32 

- 

±0. 22x10"® 

±0. 29x10"^ 


(b) Parameters for best-fitting pulsed curve 


AU series 

y-axis zero offset, 

Total particle 

Kinetic 

Floating potential, 

Effective escape cone 

Mean square error, 

run number 


flux reaching 

temperature, 

Vp 

angle in velocity space, 

percent 


A 

analyzer. 

Vo. 

V 

^o. 




ii. 

eV/ion 


deg 




A 





1 

-1. 40x10“® 

43. 7x10"® 

760 

-118 

37 

1. 43 

3 

-1. 19 

55. 2 

785 

-91. 5 

53 

2. 37 

5 

-. 69 

48. 2 

694 

-27.4 

44 

. 85 

7 

-.89 

52. 3 

730 

-85. 9 

50 

.73 

9 

-1. 30 

47.1 

792 

-196. 3 

43 

1. 14 

11 

-1.76 

49. 5 

835 

-180 

45 

1. 53 

13 

-.87 

43.8 

655 

31.2 

37 

1.37 

15 

-1.02 

44.8 

703 

-67. 5 

39 

.83 

17 

-.96 

46.2 

736 

-95.9 

41 

.73 

19 

-.96 

51.4 

752 

-49.6 

48 

. 92 

21 

-1. 70 

40.7 

783 

-147 

30 

.83 

22 

-1. 82 

41. 9 

813 

-206 

33 

1.20 

23 

-1, 21 

43.3 

756 

-137 

36 

1. 03 

Average 

-1. 21x10"® 

46. 8X10"® 

753 

j 

-105 

41 

1.15 

Mean error 

±0. 29x10"® 

1 

±3. 5X10"® 

±39 

±54 

±5 

±0. 33 
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TABLE V. - Concluded, ANALYSIS OF AU RUN SERIES 


(c) Parameters for best- fitting square-wave curve 


(d) Parameters for best-fitting delta-function curve 


AU series 

y-axis zero 

Total particle 

Kinetic 

Floating 

Mean square 

AU series 

y-axis zero 

Total particle 

Kinetic 

Mean square 

run number 

offset, 

flux reaching 

temper- 

potential. 

error, 

run number 

offset, 

flux reaching 

temper- 

error, 



analyzer, 

ature. 

vi. 

percent 


lo. 

analyzer, 

ature, 

percent 


A 

h’ 


V 



A 






A 

eV/ion 





A 

eV/ion 


1 

1. 77x10"® 

41, 6x10“® 

873 

-275 

3. 86 

1 

4. 32x10"® 

22, 9x10"® 

800 

13.1 

3 

1. 94 

43. 5 

647 

8.8 

4.93 

3 

4.73 

28.9 

700 

13.8 

5 

2. 16 

41.6 

655 

7. 29 

4.10 

5 

5. 23 

28.4 

650 

13.1 

7 

2. 08 

41.7 

624 

9.3 

4. 32 

7 

4. 92 

28.0 

650 

12.9 

9 

1. 87 

41.4 

850 

-290 

4. 49 

9 

4. 47 

22.4 

750 

13.3 

11 

1. 57 

42.1 

832 

-209 

4.12 

11 

4. 17 

23.9 

800 

12.9 

13 

2. 06 

41.6 

708 

-45 

3.65 

13 

4.61 

26.7 

750 

12.8 

15 

1.76 

41.9 

783 

-165 

3. 99 

15 

4. 39 

24. 6 

750 

13.3 

17 

1.94 

41.6 

779 

-154 

4.14 

17 

4. 59 

24.6 

750 

13.2 

19 

2. 30 

43.0 

654 

13.9 

4. 00 

19 

5. 16 

28.8 

700 

12.8 

21 

1.77 

41.4 

1016 

-464 

3. 98 

21 

4. 59 

20.9 

800 

13.7 

22 

1. 54 

41.7 

1046 

-526 

3. 99 

22 

4. 30 

20. 5 

800 

13.6 

23 

1. 93 

41. 3 

1 888 

-307 

4. 00 

23 

4. 41 

22.2 

800 

13. 2 

Average 

1. 90x10"® 

41. 9X10"® 

797 

-184 

4. 12 

Average 

4. 61X10"® 

24. 8x10"® 

746 

13.2 

Mean error 

±0. 17x10"® 

±0. 5x10"® 

±112 

±148 

±0. 21 

Mean error 

±0. 25x10"® 

±2. 6x10"® 

±44 

±0.3 


integrated energy distributions considered. As may be seen, the pulsed curve of equa- 
tion (43) gave a best fit (lowest mean square error) to the experimental data. 

The characteristics of the best-fitting curves of each type (steady, pulsed, square 
wave, and delta function) for each of the four series of runs are listed in tables HI to VI. 
In these tables the best-fitting values of I^, Ij, V^, and 0^ are listed for each in- 
dividual run. The mean square error for a given run is listed, as are the elapsed time 
and the vacuum tank pressure appropriate to each run. By examining the run-to-run var- 
iations of the parameters in tables IH to VI, one may obtain a qualitative feeling for the 
precision with which V^, 6 ^, etc., can be evaluated by the curve-fitting method used in 
this report. The greatest run-to-run variation occurs in the parameter I^. This param- 
eter, however, is not a significant experimental variable, since it merely measures the 
precision with which one can zero the dc microvolt-ammeter and/or the X-Y recorder. 
The zero offset represented by 1^ was generally about 1 percent of the collector cur- 
rent Ij^. 

The AS and AT series of runs listed in tables HI and IV both were run with helium 
gas, both had a very small variation of background pressure during the run series, and 
both gave a best fit for the pulsed production of particles. The AS series is notable in 
having a very small run-to-run variation of the best-fitting parameters, while the AT 
series has the lowest observed mean square error (and hence the best fit), for each in- 
dividual run. 
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TABLE VI - ANALYSIS OF AV RUN SERIES 


(a) Parameters for best-fitting steady curve 


AV series 

y-axis zero 

Total particle 

Kinetic 

Floating 

Effective 

Mean square 

Elapsed time. 

Tank pressure 

run number 

offset, 

Io> 

A 

flux reaching 
analyzer, 

h> 

A 

temper- 

ature, 

V 

eV/ion 

potential, 

'^1- 

V 

escape cone 
angle in 
velocity 
space, 

»o. 

deg 

error, 

percent 

min 

torr 

N/m^ 

1 

6. 53X10"® 

181X10"® 

196 

126 

29 

1.68 

0 

4.6X10"® 

6, 12X10"“^ 

3 

7. 57 

181 

193 

130 

30 

1.77 

3 

4.6 

6. 12 

5 

6. 59 

181 

196 

129 

29 

1.78 

7 

4.75 

6.22 

7 

6. 62 

185 

200 

136 

31 

1.74 

10 

4. 90 

6. 52 

9 

5. 91 

185 

206 

136 

31 

1.79 

14 

4.95 

6. 58 

11 

5. 84 

174 

181 

136 

25 

1.79 

17 

3.9 

5. 20 

13 

5. 98 

174 

180 

145 

25 

1.74 

21 

4.0 

5. 32 

15 

6. 77 

190 

208 

152 

33 

1. 77 

25 

5.1 

6.80 

17 

6.36 

180 

198 

140 

29 

1. 66 

28 

4. 8 

6. 40 

19 

5. 23 

178 

181 

142 

28 

1.60 

31 

3. 8 

5. 05 

21 

5. 09 

167 

170 

135 

20 

1.63 

34 

3.6 

4. 80 

23 

5. 84 

169 

171 

132 

22 

1.74 

37 

3.7 

4. 92 

25 

6. 99 

176 

186 

135 

26 

1. 57 

41 

4.6 

6. 12 

Average 

6. 26X10"® 

179X10“® 

190 

136 

28 

1.71 

- 

4. 41x10"® 

5.85x10"^ 

Mean error 

±0. 52x10"® 

±5xl0“® 

±11 

±5 

±3 

; 

±0. 07 

- 

±0. 47x10"® 

±0.625x10“^ 


(b) Parameters for best- fitting pulsed curve 


AV series 

y-axis zero offset. 

Total particle 

Kinetic 

Floating potential, 

Effective escape cone 

Mean square error, 

run number 

V 

flux reaching 

temperature, 

Vl- 

angle in velocity space. 

percent 


A 

analyzer. 

Vo. 

V 

®o» 




h> 

eV/ion 


deg 




A 





1 

6. 16x10"® 

180x10"® 

235 

169 

29 

1.81 

3 

7.24 

180 

231 

173 

29 

1. 89 

5 

6.26 

180 

235 

173 

29 

1.90 

7 

6. 22 

184 

241 

177 

31 

1. 89 

9 

5. 50 

184 

249 

179 

31 

1.91 

11 

5. 48 

173 

218 

175 

24 

1. 96 

13 

5. 60 

173 

218 

182 

25 

1. 88 

15 

6. 32 

189 

251 

193 

32 

1.98 

17 

5. 94 

180 

239 

181 

28 

1. 82 

19 

4.85 

178 

219 

179 

27 

1.78 

21 

4. 80 

166 

203 

174 

20 

1.76 

23 

5. 56 

168 

204 

173 

21 

1.88 

25 

6.63 

175 

224 

176 

26 

1.74 

Average 

5. 89X10"® 

178X10"® 

228 

177 

27 

1.86 

Mean error 

±0. 55x10"® 

±5X10"® 

±13 

±4 

±3 

±0.06 
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TABLE VL - Concluded. ANALYSIS OF AV RUN SERIES 

(c) Parameters for best-fitting square-wave curve (d) Parameters for best-fitting delta- function curve 


AV series 

y-axis zero 

Total particle 

Kinetic 

Floating 

Mean square 

AV series 

y-axis zero 

Total particle 

Kinetic 

Mean square 

run number 

offset, 

flux reaching 

temper- 

potential, 

error, 

run number 

offset, 

flux reaching 

temper- 

error. 



analyzer, 

ature, 

vi> 

percent 


lo. 

analyzer. 

ature. 

percent 


A 

h’ 


V 



A 

h’ 





A 

eV/ion 





A 

eV/ion 


1 

8. 89x10"® 

185x10"® 

335 

90 

3. 09 

1 

14. 9x10"® 

146x10"® 

400 

11.1 

3 

9. 89 

186 

335 

88 

3.11 

3 

15.8 

147 

400 

11.1 

5 

9. 05 

185 

337 

90 

3.10 

5 

13.6 

141 

450 

10.2 

7 

9. 22 

189 

336 

95 

3. 04 

7 

14.0 

145 



10.4 

9 

8. 83 

185 

326 

120 

3. 20 

9 

13.9 

146 



10.6 

11 

7.85 

197 

387 

-81 

3. 30 

11 

11.9 

136 



10. 1 

13 

7. 97 

179 

332 

94 

3. 07 

13 

13.7 

143 

400 

11. 1 

15 

9.79 

192 

328 

126 

2. 94 

15 

15. 4 

152 

450 

10.8 

17 

8. 94 

185 

338 

99 

3. 01 

17 

13.9 

143 

450 

10. 5 

19 

7. 55 

180 

302 

122 

2. 84 

19 

12. 9 

146 

400 

10. 9 

21 

7. 09 

201 

395 

-52 

3.27 

21 

11.8 

135 

400 

10.6 

23 

7. 87 

200 

388 

-39 

3. 24 

23 

12.6 

137 

400 

10. 5 

25 

9.13 

180 

333 

92 

2.92 

25 

13. 5 

138 

450 

10.1 

Average 

8. 62X10"® 

188X10"® 

344 

65 

3. 09 

Average 

13. 7X10"® 

143x10"® 

427 

10.6 

Mean error 

±0. 73X10"® 

±6X10"® 

±21 

±57 

±0.11 

Mean error 

±0. 9X10"® 

±4X10"® 

±25 

±0.3 


Run series AU and AV are listed in tables V and VI, respectively. These series are 
both notable in having a rather large mean square error and a rather large run-to-run 
variation. It is shown later that this imprecision is probably attributable to the rather 
large uncontrollable variation in background pressure that occurred during these two se- 
ries of runs. Such variations in vacuum tank pressure occurred whenever neon gas was 
used, as it was during these runs. 

In figure 9 are shown four runs from the AS series, each of which has the raw data 
and one of the four types of curve fits drawn on it. An examination of these four typical 
runs reveals that the best-fitting curves are not necessarily good fits, and that the 
square-wave and delta- function curves were very poor fits. In no case during this set of 
46 runs was either the delta-function or the square-wave distribution the best-fitting 
curve. It is permissible to conclude from this that the particles from this discharge are 
much closer to a Maxwellian distribution than they are to a delta- function or square-wave 
energy distribution function. 

In figure 10 are shown one typical run from each of the four series, with the type of 
curve fit that gave the lowest mean square error for that series. As may be seen, the 
curve fits were quite good in all cases. Only figure 10(d), which had a mean square error 
of 1. 68 percent, had a really noticeable systematic departure from the raw data. This 
departure is characteristic of the runs in this series. 
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(a) steady curve. Run AS -6; y-axis zero offset, Iq, -1.53x10'^ amperes-, total particle flux reaching analyzer, 218x10”^ amperes; kinetic temperature, V 
501 electron volts per ion; floating potential, V^, -303 volts; effective escape cone angle in velocity space, 0^, 69“; mean square error. 1.49 percent. 


L50xl0'^ 



0 200 400 600 800 1000 1200 1400 1600 1800 2000 2200 2400 2600 2800 


Retarding potential, V, V 

(b) Pulsed curve. Run AS -15; y-axis zero offset, Iq, -2. 58x10"® amperes; total particle flux reaching ana iyzer, 297x10"® amperes; kinetic temperature, V, 
709 electron volts per ion; floating potential, V^, -291 volts; effective escape cone angle In veixity space, 0 q, 78°; mean square error, 0.98 percent. 

Figure 9. - Best-fitting curves of each type. 



Collector current, I, A 


150x10“^ 



(d) Delta function curve. Run AS-24; y-axis zero offset, I^, 9.7x10“^ amperes; total particle flux reaching analyzer, Ij, 74. 1x10“^ amperes; kinetic temperature, Vq, 
500 electron volts per Ion; mean square error, 12. 6 percent. 


Rgure 9. - Concluded. 


30 



Collector current, I, A 


150x10"® 




r Raw data 
/ 


I 1 
2400 4600 


_J_ 

2800 


3000 
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(b) Series AT. Run AT-1; pulsed fit; y-axis zero offset, I., -18x10*® amperes; total particle flux reaching analyzer, I,, 792x10'® a mperesj 
kinetic temperature, Vg, 464 electron volts per Ion; floating potential, Vj, 33 volts; effective escape cone angle In velocity space, 0^, 41 ; 
mean square error, 0, 53 percent. 


Figure 10. - Best-fitting curves from each series. 
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Collector current, 


-Raw data 



(c) Series AU. Run AU-15; steady fit; y-axis zero offset, 1^, -0.83x10“^ ampere; total particle flux reaching analyzer, Ii, 48x10“^ amperes; 
kinetic temperature, Vq, 620 electron volts per ion; floating potential, Vj, -237 volts; effective escape cone angle in velocity space, 0^, 45°; 
mean square error, 0.79 percent. 


-Best-fitting curve 


I 1.1 1 1 1 —t f- 1 I I 

1000 1200 1400 1600 1800 2000 2200 2400 2600 2800 3000 3200 

Retarding potential, V, V 


(d) Series AV. Run AV-1; steady fit; y-axis zero offset, 6. 5x10“® amperes; total particle flux reaching analyzer, I^, 181x10“® amperes; 
kinetic temperature, Vg, 196 electron volts per ion; floating potential, V^, 126voltS; effective escape cone angle in velocity space, Oq, 29°; 
mean square error, 1. m percent 

Figure 10. - Concluded. 
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Sources of Error for Individual Runs 


It is of interest to examine the possible sources of random and systematic error for 
an individual run and thereby place limits on the possible accuracy of this method of data 
analysis. The illustrative data already presented can be used as a framework for a dis- 
cussion of various sources of error, as limitations on the precision of the calculated 
values, and of limitations on the physical models used in the interpretation of the data. 
Even the best-fitting curves exhibit some mean square error. This can result from three 
general sources: (1) the laboratory situation does not conform in all details to the theo- 
retical model, (2) the measurements themselves have accuracy limitations and are not 
simultaneous, and (3) the sampling of the curves for computer processing is not com- 
pletely representative. 

Errors of the first kind include the following effects: 

(la) The distribution function of the ions in the discharge might not be Maxwellian 
along a radius within the escape cone in velocity space. Although one would expect this 
distribution from second law considerations (it is the most probable steady-state distri- 
bution), it is not a priori obvious that the ions in the discharge are confined long enough 
to equilibrate and become Maxwellianized. 

(lb) The ion velocity distribution may not be isotropically distributed within the es- 
cape cone, or there might not be a sharp cutoff of number density at the effective escape 
cone angle. This situation can arise in the presence of nonadiabatic effects and an elec- 
tric field between the plasma source and energy analyzer. This source of error may have 
contributed significantly to the lack of accuracy in determining the escape cone angle, and 
hence the term "effective” escape cone angle has been introduced. 

(l c) The anal 3 Ttical curves with which the raw data curves are compared assume that 
the output of particles is either sharply pulsed or constant in time. A pulsed output with 
pulses of broad half width, or a rapid succession of smeared out pulses would result in an 
intermediate type of integrated energy spectrum and complicate the interpretation of the 
data. 

(l d) In the present experiment, it was not expedient to aline the axis of the energy 
analyzer along the local magnetic field lines. These lines were at an angle of about 30° 
with respect to each other under the experimental conditions, and this certainly intro- 
duced a systematic error into the calculation of the escape cone angle 9 ^. It also would 
have made the velocity distribution asymmetrical about the analyzer axis, and this may 
have had a small effect on Ij, V^, and the other parameters. Such an effect, however, 
should have been small by comparison with the other systematic and random errors pres- 
ent in this experiment. If a more accurate value of 9 ^ were desired, it would suffice to 
design the apparatus so that the analyzer axis is alined with the magnetic field lines. 

Such alinement was deemed not essential in the present series of experiments, since it 
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was desired only to show that the computer program could repeatably converge to the 
same value of 9 ^, and that the influence on the curve shape of a particular 9 ^ could be 
reliably distinguished from small random run-to-run variations in the curve shape. An 
examination of tables in to VI reveals that the computer program did, in fact, reliably 
converge to a given escape cone angle, with a rather small mean error. 

The magnitude of the curve-fitting error associated with the preceding four factors, 
individually or collectively, is difficult to estimate. Without calculating the retarding po- 
tential curve to be expected for individual anisotropic, non-Maxwellian models, there is 
no basis for comparison with the models used. The four models used as a basis for the 
computer program predict retarding potential curves of four distinct types. However, the 
reverse is not true - a curve of a given type might be obtained from some model that was 
not considered, or was intermediate among those considered. 

Systematic errors introduced by the experimental apparatus include the following: 

(2a) The discharge conditions may vary in time during an experimental run. This is 
a potentially serious source of error with this apparatus, since an individual X-Y plot re- 
quired from 1 to 2 minutes, and the discharge conditions could not be kept absolutely con- 
stant over this duration. An inspection of figures 9(b), 10(a), 10(c), and 10(d) shows that 
the raw data curve performs small undulations about the best-fitting computer- derived 
curve. A retrace of the raw data (two successive runs traced on the same piece of graph 
paper) showed that at least part of this undulation can be attributed to time variation of 
the discharge characteristics while the retarding potential curve was being swept out. 

(2b) The problems arising from analyzer design, including Debye sheaths, secondary 
electron emission, focusing effects, etc. , have been treated extensively in the literature 
cited and are not discussed herein. Precautions were taken in the design of the present 
analyzer experiment to avoid the most obvious of such difficulties, and it is felt (based on 
the very low mean square error achieved in the series AT runs) that errors introduced by 
the analyzer design did not contribute significantly to the discrepancy between the analyti- 
cal curves and the experimental data. 

(2c) A small error may have been introduced into the data by time drift of the elec- 
tronics during a run. All the apparatus was turned on to warm up at least 3 hours before 
a series of runs was taken. The variation during a given run was observed to occur only 
in the zero setting of the Y-axis of the recorder (i. e. , y, and this was no more than 
about 0. 5 percent of Ij. 

(2d) The accuracy with which can be determined in a given run depends on the 
accuracy with which the X-axis of the recorder was calibrated. In the present series of 
experiments, the retarding potential was measured with a 1-percent voltmeter, and the 
X-axis checked against this standard. The X-axis deflection agreed with the meter read- 
ing to within 1 percent. Table IV(b) shows that all the runs of the pulsed fit of the AT se- 
ries had a mean square error of about 0. 5 percent. This was the best fit of the four run 
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series studied. Since an additional inaccuracy of about 1 percent may have been intro- 
duced by the X-axis calibration, it is probably fair to conclude that the best possible ac- 
curacy with which the effective ion kinetic temperature can be determined is about 
1. 5 percent with the equipment used in the present experiment. The mean square errors 
obtained during the test runs could in some cases all be assigned to these relatively minor 
factors. 

In addition to these errors, the computer processing involves the following additional 
factors: 

(3a) Some error could be introduced into the data reduction procedure (during the 
process of reading off the X-Y plots and listing evenly spaced data points in tabular form. 
Two different individuals could read the value of collector current (for a particular value 
of the retarding potential V) and obtain values that differed by no more than about 
0. 5 percent of the value of Ij. The errors introduced in tabulating the data and in key- 
punching it for use in the computer program are generally so gross that they are readily 
detected through their production of an anomalously high mean square error in the curve 
fit. 

(3b) The computer program has been tested and refined through about 1500 curve fit- 
tings to the point where it will converge to a best-fitting curve for any integrated energy 
spectrum of monotonically decreasing form in more than 99 cases out of 100. No tend- 
ency has been observed of the program to exhibit bistable behavior, in the sense of con- 
verging to significantly different values of Ij, V^, Vj, 6 ^, and for two closely simi- 
lar raw data curves. As one may see from inspection of tables III to VI, only small dif- 
ferences in these parameters result from the small run-to-run variations that occur dur- 
ing a series of runs. It is always necessary, of course, to keep in mind the fact that a 
best fit to the experimental curve is not necessarily a good fit. 


Sources of Error for a Series of Runs 

One may obtain a useful idea of the precision of this method and of the run-to-run 
variations occurring under typical laboratory conditions by inspecting tables ni to vn. 
The second last row of tables m to VT is the numerical average of the best-fitting param- 
eters in that column, and the last row contains the mean error of this parameter, calcu- 
lated by averaging the absolute values of the differences between the average value and 
the values for the individual runs. A comparison of the mean error with the average 
value may be used as a means of estimating the run-to-run variation. In table vn is 
shown the ratio of the mean error to the average value for the four run series considered 
and for the four theoretical distribution functions considered. 

An examination of table vn shows that during run series AS and AT, the background 
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TABLE m - SUMMARY OF RUN-TO-RUN PERCENTAGE ERRORS 


Run 

series 

N 

N 

N 

Ei^vj 

N 

N 

Ei“oi 

N 

2 l^errorl 

N 

Z 

Remarks 


7 


7 

7 

7 

7 

7 



Nlo 

NIo 

NVo 

NVj 

N0O 

error 

pip 



percent 

percent 

percent 

percent 

percent 

percent 

percent 






Steady 





AS 

60 

7 . 5 

1.4 

5.6 

6.7 

7 

0 . 9 


AT 

20 

2.3 

5.9 

18 

4.8 

7.2 

.4 


AU 

30 

8.0 

6.0 

32 

13 

29 

12 

Best fit 

AV 

8.3 

2.8 

5 . 8 

3.7 

11 

4.1 

11 

Best fit 

Pulsed 

AS 

21 . 5 

13 . 9 

1.4 

9.2 

7 . 5 

11.4 

— 

Best fit 

AT 

15 

2. 1 

5.2 

15 

2 . 5 

15 

— 

Best fit 

AU 

24 

7 . 5 

5.2 

52 

12 

29 

— 


AV 

9 . 4 

2 . 8 

5.7 

2 . 3 

11 

3.2 

— 


Square wave 

AS 

7.35 

0 . 9 

3 

18.3 

— 

1.25 

- — 


AT 

5.2 

.14 

1.7 

48 

— 

1.1 

— 


AU 

9 

1.2 

14 

80 

— 

5 . 1 

— 


AV 

8 . 5 

3.2 

6.1 

88 

— 

3.6 

— 


Delta function 

AS 

4.35 

2.3 

1.8 

— 

— 

0.8 

— 


AT 

5 . 5 

1.4 

2 

— 

— 

.7 

— 


AU 

5.4 

11 

5 . 9 

— 

— 

2.3 

— 


AV 

6.6 

2.8 

5.9 

— 

— 

2.8 

— 



pressure varied by less than 1 percent of its mean value from run to run, while during 
run series AU and AV, the background pressure varied by about 12 percent. It is of 
some interest to see whether the rather large variation of background pressure might 
have been partly or mostly responsible for the run-to-run variations of V^, Vj, Ij, 6 ^, 
and Iq. In figures 11 and 12 are shown the parameters Ij and plotted as functions 
of the background vacuum tank pressure for a steady curve fit to the series AV runs. It 
is clear in this case that increasing values of and 1^^ are correlated with larger 
background pressures, although significant run-to-run variations still remain when the 
effects of changing pressure are accounted for. Figures 13 and 14 are similar plots for 
the AU run series, in which increasing values of vacuum tank pressure are associated 
with increasing values of the total collector current Ij, but decreasing values of plasma 
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190x10'® O 



Tank pressure, pj, torr 


4.8 5.2 5.6 6.0 6.4 6.8x10'^ 

Tank pressure, py, N/m^ 

Figure 11. - Total particle flux reaching analyzer as function of 
vacuum tank pressure. Series AV; steady fit; expanded scale. 



Tank pressure, py, torr 


I 1 I I I I . 

4.8 5.2 5.6 ^‘0 9 6.8x10’^ 

Tank pressure, py, Hlrvr 

Figure 12. - Kinetic temperature as function of vacuum tank 
pressure. Series AV; steady fit; expanded scale. 


37 



740 


O 


60x10'® 

56- 


£ 52 

i— 

3 

o 

j 481 

O 

o 

■S 44 




/ 




/ 

/ 

O / 

/ 

/ o 




fe° 


O/ 


2L 

e 


700 

660 

620 

580 

540 


'"nO 

o ^ 


\ 


\ 


N 


N 


\ 


N 


40l 1 L I I II 

1.4 1.6 1.8 2.0 2.2 2.4 2.6x10'® 

Tank pressure, pj-, torr 


500 L- 
1.4 


J I I I 

1.6 1.8 2.0 2.2 
Tank pressure, pp, torr 


I I 
2.4 2.6x10'^ 


1_ 1_ I I I 

1.8 2.2 2.6 3.0 3.4x10'^ 

Tank pressure, pj-, N/m^ 

Figure 13. - Total particle flux reaching analyzer as 
function of vacuum tank pressure. Series AU; 
steady fit; expanded scale. 


I 1 1 I I . 

1.8 2.2 2.6 3.0 3.4x10'4 

Tank pressure, p^, N/m^ 

Figure 14. - Kinetic temperature as function of 
vacuum tank pressure. Series AV; steady 
fit; expanded scale. 


kinetic temperature V^. As in the AV run series, some, but not all, of the run-to-run 
variations can be ascribed to the uncontrollable changes of vacuum tank pressure in this 
apparatus. 

An examination of table vn shows that the run-to-run variation in V^, Vj, and 1^ 
of the AT run series is as great as that of the AU and AV series, even though the back- 
ground pressure variation was much smaller than that of the latter two runs. An exami- 
nation of table rv shows no systematic variation of or V^ with the small variations 
of background pressure that did occur. One must look elsewhere for an explanation of the 
variations of the AT series. In figure 15 the parameters V^ and Ij are plotted as func- 
tions of the elapsed time at which the individual runs of the AT series were taken. There 
is an obvious similarity of the curves, which justifies the conclusion that much of the run- 
to-run variation in this series was the result of a time- dependent ’’hidden variable, ” 
which was some unmeasured characteristic of the discharge or apparatus other than the 
background pressure. 

The run-to-run variation of the AS series apparently represents the best possible op- 
erating conditions with the lowest systematic error. The best-fitting parameters do not 
show any systematic trends with the small fluctuations in background pressure. Fig- 
ures 16 and 17 show that the parameters V^ and Ij exhibit no systematic variation in 
time. Hidden variables are apparently not responsible for any time-dependent systematic 
errors, and therefore it is probably justified to conclude that the systematic errors in the 


38 




Figure 15. - Total particle flux reaching analyzer and kinetic tem- 
perature as functions of elapsed time. Series AT; pulsed fit. 



Figure 16. - Total particle flux reaching analyzer as function of 
elapsed time. Series AS; pulsed fit; expanded scale. 
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Figure 17. - Kinetic temperature as function of elapsed time. 
Series AS; pulsed fit; expanded scale. 


AS series were as small as possible with this apparatus. If it is assumed that the AS run 
series had no systematic variation whatever attributable to the apparatus and that all the 
run-to-run variation for this series could be ascribed to imprecision of the data handling 
method described in this report, an upper limit can then be placed on the precision of this 
method. 

The pulsed fit to the data of the AS run series (which was the best-fitting curve type) 
shows that the effective ion kinetic temperature had a run-to-run variation of only 
1. 4 percent, which compares favorably with other methods of determining this quantity. 
The floating potential Vj was determined with a run-to-run variation of 9 percent, and 
the effective escape cone angle 6 ^ varied by 7. 5 percent. The parameter Ij is propor- 
tional to the plasma density in the discharge and varied by 14 percent during the AS run 
series. This variation was larger than that of the other three run series. Even this 
large variation, however, compares favorably with the precision of Langmuir probe 
measurements of particle density made under similar conditions. 

If the measurements described herein are typical of laboratory measurements in gen- 
eral, an examination of table VII suggests that it is reasonable to expect this method to 
yield a run-to-run variation of about 5 percent of the particle energy; 20 percent of the 
floating potential; 6 percent of Ij^, and hence relative measurements of the charged par- 
ticle density; and 10 percent of the value of the effective escape cone angle. 
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Discriminatory Ability of Data Analysis Technique 


The reader has no doubt observed that run series AS and AT yielded a better fit to 
the pulsed integrated energy spectrum than to the integrated energy spectrum appropriate 
to a steady distribution of Maxwellian particles at the analyzer. An examination of ta- 
bles in and rv shows that the pulsed fit was the better one in every run of these two se- 
ries, and that the pulsed mean square error was substantially smaller than that of the 
steady curve fit in each run. This state of affairs appears very surprising in a steady- 
state experimental situation and requires comment. 

An oscilloscope trace of the collector current is displayed in figure 18 as a function 
of time for three typical but distinct operating conditions. The time base in all three 
cases is 20 microseconds per centimeter, and the zero line is the base line between 
pulses. Figure 18(a) was taken under discharge operating conditions similar to those of 
run series AS and AT, and the ions are indeed reaching the collector in pulses. It is a 
peculiar characteristic of this discharge that particles are lost from it in periodic bursts 
under a wide variety of operating conditions. The pulses occur with a periodic frequency 
that has been analyzed and reported elsewhere (refs. 22 and 23). 

Figure 18(b) shows the waveform of the current to the collector under intermediate 
operating conditions in which the pulse width is broader than the time of flight across the 
discharge, and the pulses are closer together in time. Figure 18(c) shows the collector 
current waveform under discharge operating conditions typical of run series AU and AV. 
The zero level is defined by the lowest excursion of the signal in this case. This wave- 
form is above its zero level most of the time, but it does have individual peaks above the 
constant background current. The current waveform therefore has characteristics of 
both pulsed and steady- state operation. The steady and pulsed integrated energy spectra 
are about equally good fits to the raw data in run series AU and AV, as may be seen by 
examining the mean square errors for the individual runs of these series in tables V 
and VI. In contrast, the pulsed curve fit was unquestionably the better fit in series AS 
and AT, where the current waveform was clearly pulsed in nature. 

The ability of the computer program to converge on the pulsed distribution, in a sit- 
uation in which one might have expected the steady-state distribution to be appropriate, 
is an encouraging indication that it is possible to use the raw data to discriminate cor- 
rectly among the possible energy distribution functions. 

Large changes in the effective escape cone angle will result in only relatively small 
alterations in the geometric shape of the integrated energy spectrum. It is therefore 
proper to ask whether systematic variations in the escape cone angle can be discrimi- 
nated from noise, run-to-run variations, or changes in the shape of the energy distribu- 
tion function due to changes in the other discharge characteristics. Two series of runs 
(designated CK-CL and CM-CN) were made to indicate (but not prove) that such discrim- 
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(c) High densities; low ion energies. 
Figure 18. - Collector current as function of time. 
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Figure 19. - Effective escape cone angle as function of 
magnetic field strength. Run series CK-CL 
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Figure 20. - Kinetic temperature as function of maximum 
magnetic field strength. Run series CK-CL 
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Figure 21. - Effective escape cone angle as function of maximum magnetic field strength. Run 
series CM-CN. 



ination in the effective escape cone angle is possible under laboratory conditions. 

The individual runs in these series were made at different magnetic field strengths, 
while all other discharge parameters shown in table II remained the same. The results 
are plotted in figures 19 to 22. Figures 19 and 21 show the best-fitting effective escape 
cone angle plotted as a function of magnetic field strength. There is clearly a systematic 
trend in the data of increasing escape cone angle with increasing magnetic field strength. 
The run-to-run variation of 9 ^ about a smooth curve through the data points is approxi- 
mately the ±10 percent variation that the four earlier run series led one to expect. The 
ion kinetic temperature is plotted as a function of magnetic field strength on figures 20 
and 22. A systematic trend is evident here also. The kinetic temperature and escape 
cone angle increase together to the point at which the kinetic temperature is a maximum, 
after which the kinetic temperature decreases and the effective escape cone angle is con- 
stant with magnetic field strength. 

An interpretation of these data in terms of processes occurring in the discharge is 
beyond the scope of this report. The increase of escape cone angle with magnetic field 
shown in figures 19 and 21 is surprising. However, qualitatively similar behavior was 
observed in the two run series CK-CL and CM-CN, which were taken under different dis- 
charge operating conditions, and this is considered to reinforce the contention that fig- 
ures 19 to 22 show the actual behavior of the parameters indicated. 
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CONCLUSIONS 


Analytical expressions have been derived for the integrated energy spectra which re- 
sult from a pulsed and steady- state flux of Maxwellianized particles. These energy spec- 
tra include as parameters the floating potential Vj of the source of the particles, the 
kinetic temperature of the particles, and the half-angle 9 ^ of the cone which they 
occupy in velocity space. A computer program was written which obtains best-fitting 
values of these parameters for a given experimental energy spectrum. The raw data 
from a series of experimental runs with a retarding potential energy analyzer were com- 
pared with these analytical expressions by this computer program, and the best-fitting 
values of V^, Vj^, and 6 ^ obtained for the experimental conditions. 

Certain series of runs were best fitted by the integrated energy distribution function 
appropriate to a pulsed output of particles. This result, rather surprising in the 
’’steady- state” discharge under investigation, was consistent with the results of addi- 
tional investigation. The dynamics of the discharge were such that particles were, in 
fact, reaching the energy analyzer in pulses. 

Four series of at least 10 individual runs were made to obtain typical values of the 
run-to-run variation to be expected by using this technique. Each of the four series was 
held under conditions as nearly constant as the experimental apparatus would allow. For 
the experiment used as an example, this data analysis technique typically yielded a run- 
to-run variation of about 5 percent of the particle kinetic temperature; 20 percent of the 
floating potential; 6 percent of Ij, and hence on relative measurements of the charged 
particle density; and 10 percent of the value of the effective escape cone angle. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, April 8, 1968, 

129-02-03-05-22. 
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APPENDIX A 


NOMENCLATURE 


Mathe- 

matical 

symbol 

A 

■o 

max 



V 


FORTRAN Description 

symbol 

o 

— area of energy probe aperture, m 

— maximum magnetic field, T (kG) 

— magnetic field existing at collector, T (kG) 

-1 Q 

— electronic charge, 1.602x10 C 

— geometry factor 

lA y-axis zero offset on energy spectra curve 

IB total current to analyzer, defined by eq. (49) 

— current, A 

— total current, defined by eq. (28) 

2 

— current density. A/m 

— Boltzmann constant, 1. 38x10"^^ J/°K 

— number of particles emitted by source over hemisphere in ve- 

locity space, per pulse 

Q 

— particle number density, particles/m 

O 

— particle number density of source, particles/m 

— kinetic temperature, °K 

VA particle kinetic temperature, eV 

VB floating potential of plasma, V 

VC potential defined in fig. 6(a), V 

— particle speed, m/sec 

— most probable speed, m/sec 

particle velocity in z -direction, m/sec 

— velocity defined by eq. (9), m/sec 

— average speed of particles, VSyRT/n-M, m/sec 
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Mathe- 

FORTRAN 

Description 

matical 

symbol 

symbol 


W 

— 

number of particle pulses per second 

^0 

— 

dimensionless parameter defined by eq. (33) 

e 

— 

angle with respect to z-axis and/or magnetic field line 


TH 

effective escape cone angle in velocity space 

<p 

— 

rotational angle in plane normal to z-axis and/or magnetic 
field line 
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APPENDIX B 


COMPUTER PROGRAM 
by Loretta R. Ellis 

A computer program was written to find values of the constants 1^, Ij, V^, Vj, 0^ 
in four analytical expressions, such that a best-fitting curve is obtained. These expres- 
sions are given in equations (40) to (48). The least squares method was used for each 
curve fit. This method states that the best-fitting curve is that for which the sum of the 
squares of the residuals is a minimum. To use this method conveniently, the equations 
should be linear in the unknown constants. The given expressions are not of this form; 
therefore, new expressions were derived which were linear in corrections for approxi- 
mate values of the constants. 

According to Scarborough (ref. 24), this can be done as follows: Let y = f(x, a, b) be 

a given function from expressions (40) to (48). Let a^ and b^ be approximate values of 

the unknown constants a and b, and let a and i3 denote corrections to a and b 

0 0 

such that a = a^ + O' and b = b^ + /3. Then y* = f(x, a^, b^) will be a function whose 
graph approximates the graph of y = f(x, a, b). If the graph of y = f(x, a, b) passes as 
near as possible to each of the n points (xj, y^), (X 2 ,y 2 ), . . (x^^, y^^) on the experi- 

mental curve, the residuals Vj^ will be the difference between the computed and ob- 
served values 


Vk = H\, a, b) - 


(Bl) 


or 


Vjj = i(\, a^ + O', b^ + /3) - yj^ (B2) 

where yj^ = y(xj^) represents the k*"" point on the experimental curve (k = 1, 2, , , . , n). 
Using Taylor’s expansion for a function of several variables, one obtains 


f(x, , a + Of, b + |3) = f(x. , a , b ) + Q!(— ) + /3(— ) + Higher order terms 

^ ° ° k 0 o \d^/ \9b/,_. 


(B3) 


a=ac 

b=b. 


a=a^ 

b=b. 


Omitting the higher order terms gives equation (Bl) in the form 
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I 


or 




= i(\, 






(B4) 



(B5) 


Let - Y^r) where is the residual for the approximation curve = f(Xj^, 

since it is the difference between the observed ordinates and the ordinates of this curve. 
Then 


Vk = « 


a=a 

c 

b=b^ 


9b, 


+ r. 


a=ac 

b=b. 


(B6) 


The method of least squares can be applied now to find the corrections a and given 
approximate values a^ and b^. 

In finding the best representative curves for the four expressions, better results 
were obtained by solving for two constants at a time and holding the remaining constants 
fixed, rather than trying to solve for more than two constants. Those constants which 
were the hardest to approximate were the first on which to be iterated. 

Initial values for the constants 1^, Ij, V^, and Vj are supplied by the user in the 

main program. The initial value of 9^ is taken as 40° (0^ is converted to radians for 
the computations). These initial values are used as approximate starting values in each 
of the four curve fits. 

The main program SPECTR (fig. 23) reads the input data, supplies this information 
to the subprograms, stores the results until all the calculations are completed for a 
given set of data, then prints those results. 

Subroutine MODE (fig. 24) finds the best-fitting values for the constants 1^, Ij^, V^, 
Vi, 0Q for the steady-state and for the pulsed-state expressions. Subroutines STEADY 
(fig. 25) and PULSED (fig. 26) supply the derivatives and residuals for those expressions 
for the least squares evaluation of the corrections needed (steps (1) and (2) following). 
Subroutine CALC (fig. 27) provides the necessary calculations for step (3). The itera- 
tion cycle consists of three steps: 
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(1) Maximum of five iterations to find and Vj, keeping I^, Ij, and 0^ fixed 

until the corrections are small enough 

(2) Maximum of five iterations to find and 6 ^, keeping 1^, V^, and Vj fixed 

until 

(a) Corrections are small enough 

(b) |a0q| >0.5 

(c) (0Q + A0^) falls outside the range (0, tt/2) 

(3) Maximum of 50 iterations to find keeping 1^^, V^, Vj, and 0^ fixed 
The values of the constants are updated after each iteration. On completion of 

step (3), the values are compared with those supplied at the beginning of the cycle. If 
they do not agree within the limits specified, the preceeding three steps are repeated for 
a maximum of 50 cycles or until the values do agree to three significant places. The 
values at the end of a cycle become starting values at the beginning of the next cycle. 

Subroutine SQUARE (fig. 28) finds the best fitting values for the constants 1^, Ij, 

Vq, and for the square wave function. The method used here is divided into two 
parts: 

(1) Maximum of five iterations to find and Vj, keeping 1^ and Ij fixed until 

the corrections are small enou^ 

(2) Maximum of five iterations to find 1^ and Ij, keeping and Vj fixed, until 

the corrections are small enough 

The values of the constants are now compared with those from the beginning of the cycle. 
If they do not agree within the limits specified, the preceding two steps are repeated for 
a maximum of 50 cycles or until the values do agree to three significant places. As in 
MODE, the values at the end of one cycle become starting values at the beginning of the 
next cycle. 

Subroutine FDELTA (fig. 29) finds the best-fitting values for the constants I^, Ij, 

and V for the delta-wave function. Values of I and I^ are found as V takes on 
o o 1 0 

successive values of V from the input curve. When the mean square error between the 
computed curve and the input curve reaches a minimum, the iteration stops. The maxi- 
mum number of iterations is the same as the number of points on the given curve. Con- 
trol is then returned to the main program for storage of output. 

Subroutine LESQ (fig. 30) finds the normal equations according to the least squares 
method, given the necessary derivatives and residuals for each data point. It then uses 
these equations to solve for the unknowns. The normal equations here are two equations 
in two unknowns, the unknowns being the corrections. 

Subroutine STEADY, PULSED, and CALC supply the necessary computations (deriv- 
atives, residuals, etc. , ) for subroutine MODE. 

The input consists of the following cards for each set of data: 

(1) The first card contains the run number RUN, the number of points N for this 
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particular curve, and a control indicator NV. RUN is read in alphameric format; 

N and NV are read in integer format and both must be less than 100, the maximum size 
of the V- and I-arrays. These three values RUN, N, and NV, must end in card col- 
umns 5, 10, and 15, respectively. 

(2) If NV s 0, the V-array is unchanged. If NV > 0, the first NV elements of the 
V-array are read. In order to minimize card input, the maximum number of elements 
needed of the V-array are read in, then left unchanged until new values are read again 
(when NV > 0). Hence, NV should be as large as the largest N which will be associated 
with it. 

(3) Next, the first N elements of the I-array are read in. 

The elements of the V- and I-arrays are read in floating-point format with 10 card 
columns allowed for each value. There are eight numbers per card. The last card of 
each array may have less than eight numbers on it. 

(4) The last card of the input set contains the initial values for I^, 1^, V^, and Vj^. 
The values are read in floating-point format with 10 card columns allowed for each value. 

The program stops when there are no more data to be read. Sample input is included 
with the FORTRAN listing. 
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I KTR = KTR + 1 


/ Call STDYl \ 
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Figure 24. - Flow chart for subroutine MODE. 























Figure 25. - Flow chart for subroutine STEADY. 



Figure 26. - Flow chart for subroutine PULSED. 



Figure 27. - Flow chart for subroutine CALC. 























Figure 29. - Flow chart for subroutine FDELTA. 


C 


LESQ 


) 




Finds normal equations 
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Figure 30. - Flow chart for subroutine LESQ. 
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SIBFTC SP6CTR LIST, DECK 

DIMENSION COEF(5) tZI ( 100, 5 ) , Z ( 10 , 5 ) , ZX ( 5 ) ,CO ( 5 ) 

COMMON /ALL/N,V( 100) ,I ( 100) , ICALC ( 100 ) , DEL < 100, 5 ) , lA, IB,VA, 

1 VB,TH,A,B,E,M,KOUNT,KTR 

EQUIVALENCE ( COEF ( 1 ) , I A ) , ( Z X ( 1 ) , Z I A ) , ( ZX ( 2 ) , Z I B ) , ( ZX ( 3 ) , ZV A ) , 
1 (ZX(4),ZVB),(ZX(5),ZTH) 

EQUIVALENCE (C0(1 ) , lASAV) , ( CO ( 2 ) , I BSAV ) , ( CO I 3 ) , VAS AV ) , 

1 (C0C4) ,VBSAV) , <C0(5) ,THSAV) 

REAL I , ICALC, lA, IB, lASAV, IBSAV 
DATA ZTH,ZZTH/40. ,0.69813/ 


RUN IS THE RUN NUMBER 

(V,I) ARE THE (X,Y) COORDINATES - VOLTAGE VS CURRENT 
NONE OF THE INITIAL GUESSES Z I A , Z I B , Z VA , Z VB , SHOULD i 
N = NUMBER OF POINTS (V,I) - N.LE.lOO 
NV = CONTROL INDICATOR 

NV.LE.O PREVIOUS VALUES OF V(J) ARE TO BE USED 


SHOULD BE EQUAL 


NV.GT.O NV ELEMENTS OF 
AND NV.LE.lOO 
V(J) VOLTAGE INPUT 

I(J) CURRENT INPUT 

ZI(J,K) CALCULATED CURRENT 
Z(L,K) 


TO BE USED.., 


.GE*MAXIMUM N 


USING FINAL VALUES OF IA,IB,VA,VB 


J = 1,N 
K^l,4 


L=l,10 


(DENOTES WHICH POINT) 

(K = l DENOTES STEADY-'STATE MODE) 

(K=2 DENOTES SQUARE-WAVE FUNCTION) 
(K=3 DENOTES DELTA FUNCTION) 

(K=4 DENOTES PULSED MODE) 

(L=l DENOTES lA (1(0)) ) 

(L = 2 DENOTES IB (1(D) ) 

(L=3 DENOTES VA (V(0) ) ) 

(L=4 DENOTES VB ( V( 1 ) ) ) 

(L=5 DENOTES THETA ) 

(L=6 DENOTES THE NUMBER OF ITERATIONS 
(L=7 DENOTES THE MEAN SQUARE ERROR ) 


WRITE (6,90) 
READ (5,100) 
IF (NV.GT.O) 
READ (5,110) 
READ (5,110) 


IF (ZIA.EQ.O.O) 
IF (ZIA.EQ.O.O) 
IF (ZIB.EQ.O.O) 
IF (ZVA.EQ.O.O) 
IF (ZVB.EO.O.O) 
FLN=N 

DO 20 J=l,5 
DO 20 K=l,10 
Z(K,J)=0.0 
DO 80 L=l,4 


RUN,N,NV 

READ (5,110) (V(J 
( I ( J) ,J = 1,N) 
ZIA,ZIB,ZVA,ZVB 
0) ZIA=I (N)^O.OOl 
0) ZIA=1.0E-3 
0) ZIB=I(1) 

0) ^VA=V(2) 

0) ZVB=V(2) 


(V( J),J=1,NV) 


SPE 90 
SPE 100 
SPE 110 
SPE 120 
lESPE 130 
SPE 140 
SPE 150 
SPE 160 
SPE 170 
SPE 180 
SPE 190 
SPE 200 
SPE 210 
SPE 220 
SPE 230 
SPE 240 
SPE 250 
SPE 260 
SPE 270 
SPE 280 
SPE 290 
SPE 300 
SPE 310 
SPE 320 
SPE 330 
SPE 340 
SPE 350 
SPE 360 
SPE 370 
SPE 380 
SPE 390 
SPE 400 
SPE 410 
SPE 420 
SPE 430 
SPE 440 
SPE 450 
SPE 460 
SPE 470 
SPE 480 
SPE 490 
SPE 500 



non 


30 


40 

50 

60 


70 


80 


90 

100 

110 

120 

130 


140 


IA=ZIA SPE 

IB=ZIB SPE 

VA=ZVA SPE 

VB=:ZVB SPE 

GO TO (30t40f 50»30) tL SPE 

TH=ZZTH SPE 

CALL MODE (L) SPE 

Z(5tL)-57.295779=««TH SPE 

GO TO 60 SPE 

CALL SQUARE SPE 

GO TO 60 SPE 

CALL FDELTA SPE 

S=0*0 SPE 

DO 70 J=1,N SPE 

ZI( JtL)=ICALC(J) SPE 

S=S+( I ( J )-ICALC( J ) )**2 SPE 

Z(ltL)=IA SPE 

Z(2,L)=IB SPE 

Z(3tL)=VA SPE 

Z(4tL)=VB SPE 

Z(6,L)=K0UNT SPE 

Z(8,L)=I8/AS0RT(VA) SPE 

Z(7,L)=100./I ( 1 )=J^SQRT(S/FLN) SPE 

WRITE (6t 120) RUN, ( J,V( J) T I(J)t(ZI(J,L)tL=lt4),J=l,N) SPE 

WRITE (6,130) (ZX(J),(Z(JtL),L=l,4),J=l,5),((Z(J,L),L=l,4),J=6,7) SPE 
WRITE (6,140) (Z(8,L),L-lt4) SPE 

GO TO 10 SPE 

SPE 

SPE 

SPE 

FORMAT (IHl) SPE 

FORMAT (A5,2I5) SPE 

FORMAT (8E10.0) SPE 

FORMAT (I6H RUN NUMBER , A5/ IHK , 1 2X , IHV , 1 5X , IH I , 1 3X , 6HSTE AOY , lOSPE 
1X,7HS0-WAVE,9X,5HDELTA, 11X,6HPULSED// (I5,F12*3,4X,5E16.5)) SPE 

FORMAT ( 1HL//27X,7HINITIAL,9X,6HSTEADY, lOX , 7HS0-WAVE, 9X , 5HDELT A , 1 ISPE 
1X,7HPULSED //13X,8HI0 = , 5E 16 . 5/ 1 3X ,8HI 1 = , 5E 16 . 5/ 1 3X , 8HV0 SPE 

2= ,5E16^5/13X,8HV1 = , 5E 16* 5/ 1 3X, 8HTHETA = , 5E 16* 5// 13X, 8HK0SPE 

BUNT = ,16X,4E16*5/13X,8H£RR0R = ,16X,4E16*5) SPE 

FORMAT (13X,8HNE = ,16X,4E16.5) SPE 

END SPE 


510 

520 

530 

540 

550 

560 

570 

580 

590 

600 

610 

620 

630 

640 

650 

660 

670 

680 

690 

700 

710 

720 

730 

740 

750 

760 

770 

780 

790 

800 

810 

820 
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SIBFTC MODE. LIST, DECK 

MOD 00 


SUBROUTINE MODE (L) 

MOD 10 


DIMENSION C0(10),C0EF(10) 

MOD 20 


COMMON /ALL/N,V( 100) , I ( 100 ) , I CALC ( 100 ) , DEL ( 100, 5 ) , lA, IB,VA, 

MOD 30 


1 VB,TH,A,B,E,M,ITH,KTR 

MOD 40 


EQUIVALENCE (C0EF(1),IA) 

MOO 50 


EQUIVALENCE ( CO ( 1 ) , I ASAV ) , ( CO ( 2 ) , I BSAV ) , ( CO ( 3 ) , V ASAV ) , 

MOD 60 


1 <C0(4),VBSAV),(C0(5),THSAV) 

MOD 70 


REAL I,ICALC,IA,IB,IASAV,IBSAV 

MOD 80 


DATA C0NV/1.0E--3/ 

MOD 90 


ITH=0 

MOD 100 

10 

ITH=ITH+1 

MOD no 


IASAV=IA 

MOD 120 


IBSAV=IB 

MOD 130 


VASAV=VA 

MOD 140 


VBSAV=VB 

MOD 150 


THSAV=TH 

MOD 160 


KTR^O 

MOD 170 

20 

KTR=KTR+1 

MOD 180 


IF (L.E0.4) GO TO 30 

MOD 190 


CALL STEADY 

MOD 200 


GO TO 40 

MOD 210 

30 

CALL PULSED 

MOD 220 

40 

CALL LESQ 

MOD 230 


IF ( (ABS(A/VA).LE.CONV).AND.{ABS(B/VB).LE.CONV) ) GO TO 50 

MOD 240 


IF (KTR.GE.5) GO TO 50 

MOD 250 


VA=ABS(VA+A) 

MOD 260 


VB=VB+B 

MOD 270 


GO TO 20 

MOO 280 

50 

KTR=0 

MOD 290 

60 

KTR=KTR+1 

MOD 300 


IF (L.E0,4) GO TO 70 

MOD 310 


CALL STDYl 

MOD 320 


GO TO 80 

MOO 330 

70 

call PLSDl 

MOD 340 

80 

CALL LESO 

MOD 350 


IF {(ABS(A/IB).LE.C0NV).AND.(ABS(B/TH),LE.C0NV)) go to 90 

MOD 360 


IF (KTR.GE.5.0R.ABS(B).GT,0,5) GO TO 90 

MOD 370 


IB=IB+A 

MOD 380 


TH=TH+B 

MOD 390 


IF (TH.GT.O.O^AND.TH.LE. 1.57079) GO TO 60 

MOD 400 


IB=IB~A 

MOD 410 


TH=TH“B 

MOD 420 

90 

X==ABS( IA) 

MOD 430 


IX = 0 

MOD 440 

100 

IF (X.LT.IO.) GO TO 110 

MOD 450 


IX=IX+1 

MOD 460 


X=X/10. 

MOD 470 


GO TO 100 

MOD 480 

no 

IF (X.GE.1.0) GO TO 120 

MOD 490 


IX=IX-1 

MOD 500 


x=x=no. 

MOD 510 


GO TO no 

MOD 520 

120 

0IV=10.*«IX 

MOO 530 


IDIV=0 

MOD 540 


INDX=0 

MOD 550 


IK=0 

MOD 560 

130 

IK=IK+1 

MOD 570 


IF ( INDX.GE.50) GO TO 170 

MOD 580 


60 




CALL CALC (LtINDX) 

MOD 

590 


IF ( IK.NE.l.AND.E.GF.SE) GO TO 140 

MOD 

600 


SE = E 

MOD 

610 


SIA=IA 

MOD 

620 


GO TO 160 

MOD 

630 

140 

IF { IK^NE.2) GO TO 150 

MOD 

640 


DIV=-DIV 

MOD 

650 


GO TO 160 

MOD 

660 

150 

IF ( IDIV.GE.2 ) GO TO 170 

MOD 

670 


IDIV=IDIV+1 

MOD 

680 


0IV=ABS(DIV/10. ) 

MOD 

690 


IK = 1 

MOD 

700 

160 

IA=SIA+DIV 

MOD 

710 


GO TO 130 

MOD 

720 

170 

IA=SIA 

MOD 

730 


IN=0 

MOD 

740 


DO 180 J=lt5 

MOD 

750 

180 

IF (ABS( (COEF( J)-CO( J) >/C0( J) ) •LE.CONV) IN=IN+1 

MOD 

760 


IF ( IN.NE.5,AND* ITH.LT.50) GO TO 10 

MOD 

770 


CALL CALC (L,0) 

MOD 

800 


RETURN 

MOD 

810 


END 

MOD 

820- 
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$IBFTC MODEl LISTtDECK 


STD 

00 


SUBROUTINE STEADY 


STO 

10 


REAL I fICALCtlAt IB 


STD 

20 


DIMENSION Y( 100) ,TX( 100) 


STD 

30 


COMMON /ALL/N,V( 100),I ( 100 ) , ICALC ( 100 ) ,DEL < 100,5 ) , I A , IB , VA t 

STD 

40 


1 V8tANG,A,BtE,M^K0UNT,KTR 


STD 

50 

C 



STO 

60 

C 

DERIVATIVES WITH RESPECT TO VA 

AND VB 

STD 

70 

C 



STD 

80 

10 

IF (KTR*GT*1) GO TO 20 


STD 

90 


SA=SIN(ANG)«»2 


STD 

100 


CA=COS(ANG)*=5'2 


STD 

110 


TA=SA/CA 


STD 

120 


ZA=IB/SA 


STD 

130 

20 

ZB=ZA/VA 


STD 

140 


DO 40 J=ltN 


STD 

150 


IF (V(J).LE.VB) GO TO 30 


STD 

160 


X=( V( J)~VB)/VA 


STO 

170 


EX=EXP(-X) 


STD 

180 


EXT=EXP(-X*TA) 


STD 

190 


TXY=ZB*EX^{ 1*0-EXT) 


STD 

200 


ICALC ( J )=IA+ZA«EX*( 1.0-CA^EXT) 


STD 

210 


DEL( J ,1 )=TXY^X 


STD 

220 


DEL( J,2)=TXY 


STD 

230 


GO TO 40 


STD 

240 

30 

ICALC( J )=IA+IB 


STD 

250 


DEL( Jtl )=0.0 


STD 

260 


DEL(Jt2)=0.0 


STD 

270 

40 

DEL(Jr3)=I ( J)-ICALC( J) 


STD 

280 


GO TO 100 


STD 

290 


ENTRY STDYl 


STO 

300 

C 



STD 

310 

C 

DERIVATIVES WITH RESPECT TO IB 

AND TH 

STD 

320 

c 



STD 

330 


IF (KTR.GT*!) GO TO 70 


STD 

340 


DO 60 J=1,N 


STD 

350 


IF (V(J)*LE.VB) GO TO 50 


STO 

360 


Y(J)=(V(J)“VB)/VA 


STD 

370 


TX( J)=EXP(-Y( J) ) 


STD 

380 


GO TO 60 


STD 

390 

50 

DEL(Jrl)=l*0 


STO 

400 


DEL( Jt2 )=0*0 


STD 

410 

60 

CONTINUE 


STD 

420 

70 

CA=COS( ANG) 


STO 

430 


CB=CA*^2 


STD 

440 


SA=SIN( ANG) 


STD 

450 


SB = SA=!^SA 


STD 

460 


SC=SA^SB 


STO 

470 


CA=0»5*CA*SC 


STD 

480 


SA=SB/CB 


STD 

490 


DO 90 J=ltN 


STO 

500 


IF (V(J).LE*VB) GO TO 80 


STD 

510 


TXY=EXP(-Y( J)^SA) 


STD 

520 


DEL( Jt 1)=TX( J )^a*0~CB^TXY)/SB 


STD 

530 


OEL( J,2 )=TX( J )*( (CB+Y( J)*SB)*TXY-CB)/CA*IB 

STD 

540 

80 

ICALC(J) = IA+IB=!=0EL(J,1 ) 


STD 

550 

90 

DEL(J,3)=I (J)-ICALC( J) 


STD 

560 

100 

RETURN 


STD 

570 


END 


STD 

580- 
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$IBFTC M0DE2 LISTtDECK 


PLS 

00 


SUBROUTINE PULSED 


PLS 

10 


REAL I tlCALCflAtlB 


PLS 

20 


DIMENSION Y( 100) ,TX( 100) 


PLS 

30 


COMMON /ALL/NrV( 100) t I ( 100 ) » I C ALC ( 100 ) , DEL ( 1 00 , 5 ) , lA, IB,VA, 

PLS 

40 


1 VBtANG,A,B,E,M,KOUNTtKTR 


PLS 

50 

c 



PLS 

60 

c 

DERIVATIVES WITH RESPECT TO VA 

AND VB 

PLS 

70 

c 



PLS 

80 

10 

IF (KTR*GT.l) GO TO 20 


PLS 

90 


CA=COS(ANG) 


PLS 

100 


CB=1.0-CA 


PLS 

no 


CA2=CA**2 


PLS 

120 


ZB=IB/CB/1*772454 


PLS 

130 

20 

ZA=ZB/VA 


PLS 

140 


DO 40 J=ltN 


PLS 

150 


IF (V(J).LE.VB) GO TO 30 


PLS 

160 


X2=( V( J)-VB)/VA 


PLS 

170 


X=S0RT(X2) 


PLS 

180 


TXY = ZA^‘(EXP(-X2 )-EXP(-X2/CA2) ) 


PLS 

190 


ICALCC J) = IA+IB^M 1.0-( ERF(X)-CA=i=ERF(X/CA) )/CB) 

PLS 

200 


DEL( J,1)=X«TXY 


PLS 

210 


DELIJ,2)=TXY/X 


PLS 

220 


GO TO 40 


PLS 

230 

30 

ICALC ( J )=IA+IB 


PLS 

240 


DEL( Jtl)=0.0 


PLS 

250 


DEL( J,2)=0^0 


PLS 

260 

^0 

DEL( J,3)=I(J)-ICALC(J) 


PLS 

270 


GO TO 100 


PLS 

280 


ENTRY PLSDl 


PLS 

290 

C 



PLS 

300 

c 

DERIVATIVES WITH RESPECT TO IB 

AND TH 

PLS 

310 

c 



PLS 

320 


IF (KTR^GT.l) GO TO 70 


PLS 

330 


DO 60 J=ltN 


PLS 

340 


IF (V(J).LE.VB) GO TO 50 


PLS 

350 


Y{ J)=SORT( ( V( J )-VB)/VA) 


PLS 

360 


TX(J)=ERF(Y(J) ) 


PLS 

370 


GO TO 60 


PLS 

380 

50 

DEL( J,1 )=1*0 


PLS 

390 


DEL( Jr2)=^0*0 


PLS 

400 

60 

CONTINUE 


PLS 

410 

70 

SA=SIN( ANG) 


PLS 

420 


CA=COS( ANG) 


PLS 

430 


CB=1*0-CA 


PLS 

440 


CC = IB^SA/CB=5'=Cc2 


PLS 

450 


CD = 1*1283792^'CB/CA 


PLS 

460 


DO 90 J=1,N 


PLS 

470 


IF (V(J).LE*VB) GO TO 80 


PLS 

480 


X=Y( J) 


PLS 

490 


XX=ERF(X/CA) 


PLS 

500 


OEL( J,1 )=1.0~(TX( J )~CA«XX)/CB 


PLS 

510 


DEU J,2)=CC*(CD>l'EXP(-(X/CA)=!‘#2)*X-XX+TX( J) ) 

PLS 

520 

80 

ICALC ( J ) = IA+IB#DEL(J,1 ) 


PLS 

530 

90 

DEL(J,3)=I (J)-ICALC(J) 


PLS 

540 

100 

RETURN 


PLS 

550 


END 


PLS 

560- 
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$IBFTC CALCl LIST, DECK 

CAL 

00 


SUBROUTINE CALC ( L» INOX) 

CAL 

10 


DIMENSION TX(IOO) 

CAL 

20 


REAL ItICALC,IA,IB,IASAV,IBSAV 

CAL 

30 


COMMON /ALL/NtV( 100) ( 100) t ICALC ( 100 ) , DEL ( 100 , 5 ) , 1 A , I B , VA t 

CAL 

40 


1 VB,ANG,A,B,EtMtKOUNT,KTR 

CAL 

50 


E = 0^0 

CAL 

60 


INDX=INDX+1 

CAL 

70 


IF ( INDX^GT.l ) GO TO 50 

CAL 

80 


MM=l+L/3 

CAL 

90 


GO TO ( 10,30) , MM 

CAL 

100 

10 

SA=SIN( ANG)*»2 

CAL 

110 


CA=COS( ANG)=^*2 

CAL 

120 


TA=SA/CA 

CAL 

130 


DO 20 J=1,N 

CAL 

140 


IF (V(J).LE,VB) GO TO 20 

CAL 

150 


X={ V(J )“VB)/VA 

CAL 

160 


TX( J)=EXP(-X)/SA^( 1.0-CA*EXP(-X^TA) ) 

CAL 

170 

20 

CONTINUE 

CAL 

180 


GO TO 50 

CAL 

190 

30 

CA=COS(ANG) 

CAL 

200 


CB=1.0-CA 

CAL 

210 


DO 40 J=^1,N 

CAL 

220 


IF {V(J).LE^VB) GO TO 40 

CAL 

230 


X=SQRT( (V( J )-VB)/VA) 

CAL 

240 


TX( J )=l,0-‘(ERF(X)*CA'i'ERF(X/CA) )/CB 

CAL 

250 

40 

CONTINUE 

CAL 

260 

50 

DO 70 J=1,N 

CAL 

270 


IF (V(J)*LE*VB) GO TO 60 

CAL 

280 


ICALC( J) = IA+IB=S'TXU) 

CAL 

290 


GO TO 70 

CAL 

300 

60 

ICALC( J )=IA+IB 

CAL 

310 

70 

E=E+( I (J)-ICALC( J) 

CAL 

320 


RETURN 

CAL 

330 


END 

CAL 

340- 


SIBFTC LESQl LIST, DECK 

LES 

00 


SUBROUTINE LESQ 

LES 

10 


COMMON /ALL/N,V( 100) , I ( 100) , I CALC ( 1 00 ) , DEL ( 100 , 5 ) , I A , I B , VA , 

LES 

20 


1 VB,ANG,A,B,E,M,KOUNT,KTR 

LES 

30 


DIMENSION SUM(5,6) 

LES 

40 


DO 20 K=l,2 

LES 

50 


DO 20 J=K,3 

LES 

60 


SUM(K, J )=0.0 

LES 

70 


DO 10 L=1,N 

LES 

80 

10 

SUMIK, J )=SUM(K, J )+DEL ( L , K ) *DEL ( L , J ) 

LES 

90 

20 

SUM(J,K)=SUM(K,J) 

LES 

100 


A = 0, 

LES 

110 


B=0. 

LES 

120 


DENOM=SUM( 1,1 )*SUM(2,2 )-SUM( 1,2)*=«=2 

LES 

130 


IF ( DENOM.EQ.O. ) GO TO 30 

LES 

140 


A=(SUM( 1,3)=^SUM(2,2)-SUM( 1,2)*SUM(2,3) )/DENOM 

LES 

150 


B=( SUM(1,1 )#SUM(2,3)-SUM( 1,3)*SUM(2,1 ) )/DENOM 

LES 

160 

30 

RETURN 

LES 

170 


END 

LES 

180- 
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$IBFTC 

SQU LISTtDECK 



SQU 

00 


SUBROUTINE SQUARE 



SQR 

10 


REAL I,ICALC,IA,IB 



SQR 

20 


COMMON /ALL/NtVCIOO) tl (1 00 ) , I C ALC ( 100 ) t DEL ( 1 00 , 5 ) , I A , I B » VA , 

SQR 

30 

1 

VBfANGtAtBrEtMtKOUNTfKTR 



SQR 

40 


DIMENSION C0EF(5) 



SQR 

50 


EQUIVALENCE (IA,COEF(D) 



SQR 

60 


DATA CONV/l.OE-3/ 



SQR 

70 


K0UNT=0 



SQR 

80 

10 

K0UNT=K0UNT+1 



SQR 

90 


SIA=IA 



SQR 

100 


SIB=IB 



SQR 

110 


SVA=VA 



SQR 

120 


SVB=VB 



SQR 

130 


DO 120 M=lt2 



SQR 

140 


KTR=0 



SQR 

150 

20 

KTR=KTR+1 



SQR 

160 


VA=ABS( VA) 



SQR 

170 


VC=2,0^VA+VB 



SQR 

180 


IF (M.EQ.2) GO TO 70 



SQR 

190 

C 




SQR 

200 

C 

DERIVATIVES WITH RESPECT TO VA 

AND 

VB. 

SQR 

210 

C 




SQR 

220 


ALPHA=IB/ VA/2.0 



SQR 

230 


KAB = 3 



SQR 

240 


DO 60 J=1»N 



SQR 

250 


IF (V(J).LE.VB) GO TO 30 



SQR 

260 


IF (V(J).GT.VC) GO TO 40 



SQR 

270 


DEL( J f2)=ALPHA 



SQR 

280 


DEU J,1 )=DEL( Jt2 )=i^( V( J )-VB)/VA 



SQR 

290 


ICALC( J )=IA+DEL( Jt2 )«(VC-V( J) ) 



SQR 

300 


GO TO 60 



SQR 

310 

30 

ICALC ( J )=IA+IB 



SQR 

320 


GO TO 50 



SQR 

330 

40 

ICALC( J)=IA 



SQR 

340 

50 

0EL(J,1)=0.0 



SQR 

350 


DEL( Jf2)=0»0 



SQR 

360 

60 

DEL(J,3)=I(J)-ICALC{J) 



SQR 

370 


GO TO 110 



SQR 

380 

C 




SQR 

390 

C 

DERIVATIVES WITH RESPECT TO lA 

AND 

IB. 

SQR 

400 

C 




SQR 

410 

70 

KAB = 1 



SQR 

420 


DO 100 J=ltN 



SQR 

430 


DEL(Jt1)=1^0 



SQR 

440 


IF (V(J) .LE,VB) go to 80 



SQR 

450 


IF (V(J).GT.VC) go to 90 



SOR 

460 


DEL( Jt2)=( VC-V( J ) )/VA/2.0 



SQR 

470 


ICALC (J )=IA+IB^DEL (J,2) 



SQR 

480 


GO TO 100 



SQR 

490 

80 

ICALC( J)=IA+IB 



SQR 

500 


DEL( J,2)=l,0 



SQR 

510 


GO TO 100 



SQR 

520 

90 

ICALC( J)=IA 



SQR 

530 


DEL( Jf2)=0*0 



SOR 

540 

100 

DEL( Jt3)=I { J)-ICALC( J) 



SOR 

550 

110 

CALL LESQ 



SQR 

560 


KAC=KAB+1 



SQR 

570 


IF ( (ABS(A/COEF(KAB) ) *LE*CONV) 

• AND. 

(ABS(B/COEF(KAC) ).LE.CONV) ) 

GO SQR 

580 
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ITO 120 

SQR 

590 


IF (KTR*GE.5) GO TO 120 

SQR 

600 


COEF(KAB)=COEF<KAB)+A 

SQR 

610 


COEF(KAC)=COEF(KAC)+B 

SQR 

620 


GO TO 20 

SQR 

630 

120 

CONTINUE 

SQR 

640 


IF ( (ABS( (SIA-IA)/SIA) .LE.CONVI.AND, (ABS( ( S I B- I B ) / S I B ) • L E .CONV ) 

•ANSQR 

650 


10*(ABS( (SVA-VA)/SVA).LE.CONV) *AND*(ABS( ( 5VB-VB ) /SVB ) .LE .CONV ) ) 

GO SQR 

660 


2T0 130 

SOR 

670 


IF (KOUNT.lt. 50) GO TO 10 

SQR 

680 

130 

RETURN 

SOR 

690 


END 

SQR 

700- 


$IBFTC FDEL LIST, DECK 

DEL 

00 


SUBROUTINE FDELTA 

DEL 

10 


REAL ItlCALCtlAflB 

DEL 

20 


COMMON /ALL/N,V( 100) , I( 100) , I CALC < 1 00 ) , DEL ( 100 , 5 ) , I A , I B t V A, 

DEL 

30 


1 VBTANG,AtB,EtMTKOUNT,KTR 

DEL 

40 


VB=0.0 

DEL 

50 


DO 10 J=1,N 

DEL 

60 


DEL(J,2)=0.0 

DEL 

70 

10 

DEL( Jtl )=1.0 

DEL 

80 


K0UNT=0 

DEL 

90 

20 

K0UNT=K0UNT+1 

DEL 

100 


VA=V(KOUNT) 

DEL 

110 


DEL(K0UNT,2)=1.0 

DEL 

120 


E=0.0 

DEL 

130 


DO 30 J=1,N 

DEL 

140 


ICALC( J )=IA+IB^DEL( Jt2) 

DEL 

150 

30 

DEL( J,3)=I ( J)-ICALC( J) 

DEL 

160 


CALL LESQ 

DEL 

170 


IA=IA+A 

DEL 

180 


IB=IB+B 

DEL 

190 


DO 40 J=1,N 

DEL 

200 


ICALC(J) = IA+IB=<=DEL(J»2) 

DEL 

210 

40 

E = E+(I(J)-ICALC(J))^=!^2 

DEL 

220 


IF (KOUNT.EQ.l ) GO TO 50 

DEL 

230 


IF (ESAVE.LE.E) GO TO 60 

DEL 

240 

50 

ESAVE=E 

DEL 

250 


SIA=IA 

DEL 

260 


SIB=IB 

DEL 

270 


SVA=VA 

DEL 

280 


IF (KOUNT.LT.N) GO TO 20 

DEL 

290 

60 

IA=SIA 

DEL 

300 


IB=SIB 

DEL 

310 


VA=SVA 

DEL 

320 


DO 70 J=1,N 

DEL 

330 


ICALC(J)=IA 

DEL 

340 

70 

IF (V(J).LT.VA) ICALC( J)=IA+IB 

DEL 

350 


RETURN 

DEL 

360 


END 

DEL 

370- 
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$OATA 

AS --3 

0 . 

57 64 

50 . 

100 . 

150 . 

400 . 

450 . 

500 . 

550 . 

800 . 

850 . 

900 . 

950 . 

1200 . 

1250 . 

1300 . 

1350 . 

1600 . 

1650 . 

1700 . 

1750 . 

2000 . 

2050 . 

2100 . 

2150 . 

2400 . 

2450 . 

2500 . 

2550 . 

2800 . 

2850 . 

2900 . 

2950 . 

141.0 

134.0 

117.8 

102.0 

63.0 

56.6 

51.75 

46.5 

27.0 

25.0 

22.5 

21.0 

13.5 

12.5 

11.0 

10.5 

7.0 

6.0 

5.5 

5.0 

3.0 

2.0 

1.0 

1.0 

^ 1.0 

- 1.0 

- 1.5 

- 1.5 

- 4.0 

- 8 . 

145 . 

350 . 

30 . 


200 . 

250 . 

300 . 

350 . 

600 . 

650 . 

700 . 

750 . 

1000 . 

1050 . 

1100 . 

1150 . 

1400 . 

1450 . 

1500 . 

1550 . 

1800 . 

1850 . 

1900 . 

1950 . 

2200 . 

2250 . 

2300 . 

2350 . 

2660 . 

2650 . 

2700 . 

2750 . 

3000 . 

3050 . 

3100 . 

3150 . 

92.3 

82.5 

75.0 

69.0 

44.25 

37.5 

33.0 

30.0 

19.0 

17.0 

16.0 

14.5 

10.0 

9.0 

8.5 

7.5 

4.5 

4.0 

3.5 

3.0 

.5 

.0 

.0 

-.5 

- 2.0 

- 3.0 

- 3.5 

- 4.0 
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